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IX 


Executive  Summary 


The  accumulation  of  space  debris  in  the  geosynchronous  region  (GEO)  has  raised  attention  among 
spacefaring  nations.  The  current  mitigation  measure  supported  and  voluntarily  practiced  is  to 
boost  satellites  into  supersynchronous  orbits  in  the  time  before  stationkeeping  fuel  is  expected 
to  be  exhausted  and  to  passivate  satellites  of  all  stored  and  generated  energy  sources  to  prevent 
explosions,  which  are  by  far  the  largest  source  of  fragmentation  debris.  Because  this  solution  does 
not  remove  mass  from  space,  debris  generation  by  other  fragmentation  events  remains  a  possibility. 
The  collision  hazard  between  inactive  satellites  in  the  supersynchronous  region  raises  questions 
about  the  consequences  of  collisions  in  this  regime  and  possible  interaction  with  GEO. 

In  considering  the  use  of  supersynchronous  orbits  for  satellite  disposal,  the  first  concern  is 
to  determine  the  minimum  safe  distance  above  GEO  such  that  objects  in  the  disposal  orbits  will 
not  interfere  with  the  GEO  population  in  the  future.  This  involves  defining  the  useful  GEO  area 
and  studying  the  perturbation  effects  on  objects  in  supersynchronous  orbits.  Thus  far,  research 
has  focused  on  propagating  the  orbits  of  intact  objects.  However,  in  the  aftermath  of  a  collision, 
pieces  of  varying  sizes  and  shapes  can  be  found  in  orbits  quite  different  from  the  parent  objects’ 
orbits. 

This  document  summarizes  background  information  on  debris  in  the  GEO  region,  sources 
and  management  strategies,  and  then  addresses  the  following  questions:  Will  orbits  of  fragments 
from  a  collision  in  a  storage  orbit  occupy  GEO  altitudes  at  some  time  after  the  collision?  If  so,  at 
what  altitude  should  the  storage  orbit  occupy  such  that  collision  fragments  will  not  interfere  with  the 
GEO  population?  The  methods  and  tools  by  which  the  effects  of  collisions  in  the  supersynchronous 
region  can  be  analyzed  are  discussed.  A  low- velocity  collision  model  is  employed  to  provide  delta- 
velocities  imparted  to  the  fragments.  An  analytical  study  of  perturbation  effects,  including  solar  and 
lunar  third  body  gravitation,  Earth  oblateness  through  degree  and  order  four,  and  solar  radiation 
pressure,  follows  in  order  to  evaluate  the  magnitude  of  these  disturbing  forces  on  the  fragmentation 
debris.  Validation  of  these  results  by  numerical  analysis  using  proven  numerical  and  semianalytical 
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orbit  propagators  is  discussed.  The  results  show  that  currently  practiced  reorbiting  distances  above 
GEO  do  not  isolate  debris  from  GEO  after  the  occurrence  of  collisions  in  the  storage  orbit. 

This  report  is  based  on  the  dissertation  work  of  the  primary  author  which  was  funded  by 
the  Air  Force  Palace  Knight  Program.  The  opinions  expressed  in  this  document  are  those  of  the 
author  and  do  not  reflect  those  of  the  Air  Force,  the  Department  of  Defense,  nor  the  United  States 
government. 
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Chapter  1 


Introduction 

Space  debris  can  be  defined  as  “any  man-made  Earth-orbiting  object  which  is  non-functional 
with  no  reasonable  expectation  of  assuming  or  resuming  its  intended  function  or  any  other  function 
for  which  it  is  or  can  be  expected  to  be  authorized,  including  fragments  and  parts  thereof.”  [29] 
Specifically,  orbital  debris  include  [46] 

•  Nonfunctional  spacecraft 

•  Rocket  bodies 

•  Mission-related  debris 

1.  Exhaust  products 

2.  Objects  released  in  spacecraft  deployment  and  operations 

3.  Refuse  from  manned  missions 

•  Fragmentation  debris 

1.  Explosion  fragments 

2.  Collision  fragments 

3.  Products  of  deterioration 

The  term  orbital  debris  is  often  used  interchangeably  with  space  debris,  as  is  the  case  in 
this  writing.  It  is  a  subject  that  is  gaining  public  recognition  in  recent  years,  with  short  articles 
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now  appearing  even  in  local  newspapers.  Despite  its  novelty  status  in  the  public  conscience,  orbital 
debris  has  been  a  topic  of  research  for  decades  now.  Space  agencies  regard  it  as  a  space  hazard  which 
previously  included  such  things  as  meteoroids  and  radiation,  problems  to  be  studied  and  to  find 
protection  against.  Unlike  the  other  space  hazards  mentioned,  however,  debris  is  man-generated, 
and  therefore,  to  a  certain  extent,  the  debris  situation  is  ours  to  assess,  control,  and  counteract. 
The  background  information  on  orbital  debris  is  extensive,  as  it  concerns  all  regions  of  near-Earth 
space  currently  in  use.  The  rest  of  this  introduction  briefly  summarizes  the  debris  situation  in  the 
geosynchronous  Earth  orbits  (GEO)  region.  This  is  followed  by  discussions  of  debris  sources  and 
proposed  management  schemes. 


1.1  Space  Debris  in  GEO 

The  focus  of  this  research  is  on  the  GEO  debris  environment,  which  differs  on  significant 
points  from  the  situation  in  low  Earth  orbits  (LEO)  as  described  in  this  section.  The  great  majority 
of  the  cataloged  debris  is  in  LEO.  Undoubtedly,  most  debris,  cataloged  as  well  as  uncataloged,  is  in 
LEO  due  to  its  extensive  use  and  the  wider  range  of  altitudes  that  comprises  LEO.  However,  at  GEO 
altitudes,  only  objects  larger  than  1  m  in  diameter  are  regularly  tracked  and  cataloged,  whereas 
10-50  cm  is  the  limit  in  LEO  [29];  therefore,  a  whole  population  of  smaller  pieces  potentially  exists 
in  GEO  undetected.  It  is  also  possible  that  fragmentation  events,  especially  if  the  parent  objects’ 
orbits  do  not  markedly  change  and  resulting  debris  are  small,  occur  unnoticed  [39]. 

Another  difference  is  that  relative  velocity  between  objects  in  GEO  averages  about 
500  m/s  [26],  compared  to  about  10  km/s  for  LEO  [49].  Since  probability  of  collision  is  a  function  of 
the  average  relative  velocity  between  objects  as  well  as  spatial  density,  greater  attention  is  directed 
at  LEO,  particularly  at  the  crowded  900-1000  km  and  1400-1500  km  altitude  bands  [46].  With 
manned  missions  currently  limited  to  LEO,  a  good  deal  of  debris  research  aims  at  evaluating  the 
risk  to  and  protecting  the  Space  Shuttle  and  the  planned  International  Space  Station  [47].  Finally, 
the  proximity  of  LEO  allows  for  active  retrieval  of  satellites  so  that  their  space-exposed  surfaces, 
which  are  records  of  impacts  from  debris,  may  be  studied.  The  result  of  these  distinctions  is  that 
there  are  more  measurement  data  and  environment  models  for  LEO  than  for  the  debris  in  GEO. 
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Then,  too,  atmospheric  drag  acts  as  a  natural  cleanser  for  parts  of  LEO.  Satellites  below 
500  km  altitude  will  deorbit  within  several  years  [32].  LEO  satellites  operating  above  500  km  can 
reduce  their  orbital  lifetime  by  reorbiting  at  end  of  life  (EOL)  to  lower  altitudes  or  simply  by 
lowering  their  perigee  as  appropriate  to  utilize  atmospheric  drag.  In  GEO,  orbital  lifetimes  exceeds 
thousands  of  years,  and  deorbiting  is  not  economically  feasible  with  current  technology.  Barring 
great  advances  in  propulsive  technology,  the  population,  intact  or  otherwise,  in  GEO  is  more  or  less 
permanent.  With  this  difference  in  mind,  debris  management  strategies  for  GEO  must  be  tailored 
to  the  situation  there.  The  rest  of  this  document  discusses  only  issues  pertaining  to  GEO. 

1.2  Physical  Characteristics  of  GEO 

The  primary  distinguishing  feature  of  GEO  is  that  satellites  placed  there  have  an  orbital 
period  of  about  1436.2  min,  one  sidereal  day.  A  distinction  can  be  made  between  the  geosynchronous 
region  and  the  particular  circular  orbits  at  about  35,786  km  altitude  and  zero  inclination,  which  are 
referred  to  as  geostationary  orbits  (GSO).  GSO  satellites  appear  to  be  fixed  in  the  sky  over  locations 
on  Earth’s  equator.  Consequently,  antennas  at  ground  relay  stations  can  be  small,  nonmoving, 
high-gain,  and  low-cost,  thus  extending  satellite  coverage  to  millions  of  direct  users  for  whom  the 
expense  of  tracking  antennas  is  out  of  reach  [61].  Satellites  in  other  GEO  orbits  have  “figure 
eight”  groundtracks,  which  grow  larger  with  higher  inclinations.  The  many  uses  of  GEO  include 
communications,  meteorology,  and  governmental  and  military  support. 

Due  to  tesseral  harmonics  effects,  satellites  must  perform  longitudinal,  i.e.,  east-west,  sta¬ 
tionkeeping  maneuvers  to  avoid  drifting  out  of  GSO.  The  ellipticity  of  the  Earth’s  equator  produces 
four  equilibrium  points  in  GSO,  two  stable  points  at  75.3°E  and  255.3°E  and  two  unstable  points 
at  165.3°E  and  345.3°E.  Uncontrolled  satellites  migrate  toward  the  stable  points.  Lunar  and  so¬ 
lar  gravity  torque  the  satellite  orbit  plane  toward  the  ecliptic  [26];  thus,  GSO  satellites  must  also 
perform  north-south  stationkeeping  maneuvers  to  control  inclination.  The  orbital  inclinations  of 
uncontrolled  satellites  increase  from  0°  to  almost  15°  in  27  years  and  then  back  to  0°  in  another 
27  years  [62].  Long  term  propagations  show  radial  stability  throughout  the  region,  defined  as 
2000  km  above  and  below  GSO  altitude  [26].  Uncontrolled,  intact  GEO  and  GSO  objects  will 
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remain  in  the  GEO  region  for  at  least  several  centuries.  Due  to  this,  the  accumulation  of  debris  is 
of  concern  for  many  GEO  operators. 

1.3  Debris  Sources 

The  first  GEO  satellite,  Syncom  3,  was  launched  in  1964.  It  relayed  signals  to  broadcast  the 
Tokyo  Olympics  to  half  of  the  Earth  [52].  Within  three  decades,  hundreds  of  satellites  have  been 
placed  into  GEO  with  an  average  of  about  25  new  payloads  per  year  at  present  [42].  Along  with 
this  population  of  active  and  inactive  payloads  is  an  assortment  of  rocket  bodies,  mission-related 
debris,  and  fragmentation  debris.  GEO  objects  can  be  divided  into  two  classes:  those  that  are  large 
enough  to  be  cataloged  and  those  that  are  not. 

1.3.1  Cataloged  Objects 

Objects  that  can  be  identified  and  regularly  tracked  are  listed  in  the  US  Space  Command 
satellite  catalog  along  with  their  two-line  element  (TLE)  sets.  This  catalog  is  not  generally  available 
to  the  public,  but  portions  of  it  can  be  obtained  through  other  sources.  From  the  NASA-issued 
TLE,  there  were  632  cataloged  objects  in  GEO  in  March  1997,  if  GEO  is  defined  by  [30] 

•  eccentricity  smaller  than  0.1, 

•  mean  motion  between  0.9  and  1.1  revolutions  per  sidereal  day,  corresponding  to  a  radius  of 
approximately  42,164  ±  2800  km,  and 

•  inclination  smaller  than  20°. 

These  objects  include  active  and  inactive  payloads,  apogee  kick  motors,  and  other  propulsive 
stages.  Additionally,  there  are  40  objects  which  are  listed  as  geostationary  or  near-stationary  in 
the  RAE  Table  of  Earth  Satellites  but  are  not  in  the  NASA  TLE  list  [30].  Objects  are  occasionally 
mislabeled  or  lost  from  tracking. 
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1.3.2  Uncataloged  Objects 

Of  the  smaller  debris  types,  mission-related  and  fragmentation  debris,  these  are  known  to 
exist  but  cannot  be  tracked.  Solid  rocket  motors  eject  large  quantities,  as  many  as  1020  in  a  single 
firing,  of  aluminum  oxide  (AI2O3)  particles.  These  particles  are  usually  no  larger  than  10  jim  in 
diameter  with  large  area-to-mass  ratios.  Thus,  they  are  strongly  affected  by  solar  radiation  pressure, 
which  causes  eventual  deorbiting  [32].  Other  mission-related  debris  too  small  to  be  detected  are 
objects  released  during  spacecraft  deployment  and  operations,  such  as  lens  covers,  shrouds,  springs, 
and  explosive  bolt  pieces.  Fragmentation  debris  also  exist.  Three  explosions  have  been  reported 
in  GEO  thus  far  [67].  One  was  a  USSR  Ekran  satellite  photographed  in  June  1978  as  it  exploded 
from  a  suspected  nickel- hydrogen  battery  failure  [39].  The  other  events  are  the  explosions  of  two 
Titan  2  transtage,  including  one  on  21  Feb  1992  that  had  been  in  orbit  since  1968  [50].  In  general, 
explosions  may  be  categorized  by  their  causes: 

•  Propulsion-related 

•  Deliberate  actions 

•  Accidental  detonations 

•  Electrical  failures 

It  is  suspected  that  more  explosions  may  have  occurred  unrecorded  or  unconfirmed  [39,  50]. 
In  LEO,  explosions  have  been  the  most  prolific  source  of  cataloged  debris.  Additionally,  long  term 
exposure  to  radiation  and  thermal  cycling  causes  loss  of  flexibility  in  plastics  and  polymers  which 
can  lead  to  structural  failures  and  gradual  fragmentation  [32].  The  large  number  of  objects  in 
GEO  and  their  long  orbital  lifetime  form  a  steady  source  pool  for  fragmentation  debris.  As  the 
number  of  satellites  and  fragments  increase,  accordingly  the  probability  of  collisions  of  every  type, 
i.e.,  satellite-satellite,  satellite-fragment,  and  fragment-fragment,  will  increase.  Some  theorize  that 
beyond  a  certain  average  debris  flux  for  a  region,  the  debris  environment  could  reach  a  state  of 
being  self-regenerative  [40],  although  this  hypothesis  is  still  the  subject  of  much  debate. 
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1.4 


Debris  Management 


There  are  no  international  treaties  or  laws  regulating  the  growth  of  orbital  debris.  However, 
various  spacefaring  groups  are  individually  drawing  up  and  enforcing  their  own  debris  mitigation 
guidelines.  A  survey  of  industry  and  civil  government  agencies  and  organizations,  conducted  by 
the  AIAA  Orbital  Debris  Study  Group  and  published  in  1992,  shows  that  several  groups  have  vol¬ 
untarily  adopted  various  methods  designed  to  mitigate  the  growth  of  orbital  debris.  This  suggests 
that  there  are  some  debris  minimizing  techniques  with  acceptable  cost-to-benefit  ratios,  an  impor¬ 
tant  criterion.  Unfortunately,  the  Study  Group  concluded  that  current  measures  are  insufficient  to 
control  the  debris  environment  [63]. 

In  the  following  sections,  active  and  passive  methods  for  orbital  debris  management  are 
described.  In  particular,  it  is  important  to  consider  the  methods’  technological  maturity  and 
economic  feasibility.  Debris  management  aims  to  minimize  risk  to  satellites  directly  and  indirectly 
through  control  of  the  growth  of  debris. 

1.4.1  Passive  Methods 

Passive  methods  are  grouped  together  by  one  common  feature:  these  methods  do  not  seek 
to  control  debris  generation.  Some  passive  methods  are 

•  shielding  of  satellites, 

•  ground-initiated  collision  avoidance,  and 

•  on-board  collision  avoidance  systems. 

Shielding  requirements  in  the  early  space-launching  days  were  based  on  the  natural  mi¬ 
crometeoroid  flux  through  near-Earth  space.  The  planned  International  Space  Station’s  shielding 
now  accounts  for  the  debris  environment  as  well  as  for  meteoroids  [47].  Likewise,  some  satellite 
designers  now  include  debris  hazard  analysis  [4].  However,  there  is  obviously  a  limit,  imposed  by 
economics  as  well  as  technology,  beyond  which  additional  shielding  is  not  feasible. 

The  second  method  uses  ground  facilities  to  propagate  the  orbit  of  a  satellite  of  interest 
as  well  as  those  of  nearby  cataloged  objects  to  predict  possible  collisions.  This  is  practiced  by 
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some  satellite  operators  [58].  Complications  arise  from  uncertainties  in  the  orbit  determination 
and  propagation  processes  as  well  as  the  satellite  operators’  desire  to  conserve  fuel  and  thus  to 
avoid  unnecessary  maneuvers.  There  is  also  some  concern  that  ground  control-initiated  avoidance 
maneuvers  cannot  be  relied  upon  since  only  about  10%  of  the  hazardous  LEO  debris  population  is 
cataloged  [65].  The  situation  is  likely  worse  in  GEO  where  typically  only  objects  larger  than  1  m 
in  diameter  are  cataloged. 

On-board  collision  avoidance  systems  are  not  technically  feasible  at  present.  There  are  two 
main  problems.  First,  tracking  objects  and  determining  their  orbits  from  a  spacecraft  present  many 
difficulties.  Secondly,  it  is  doubtful  whether  the  spacecraft  would  have  enough  time  after  detection 
of  a  possible  collision  hazard  to  maneuver  out  of  harm’s  way  [47]. 

1.4.2  Active  Methods 

Active  methods  attempt  to  directly  influence  the  debris  environment.  They  include 

•  minimization  of  space  operations  debris, 

•  minimization  of  risk  of  explosions, 

•  reorbiting  spacecraft  at  EOL,  and 

•  active  collection  or  removal. 

There  are  various  proposals  to  decrease  space  operation  debris  on  future  launches.  Simple 
solutions  include  securing  parts,  such  as  shrouds  and  lens  covers,  to  the  spacecraft  or  rocket  body. 
Other  techniques  considered  are  “bagging”  fragments  from  pyrotechnic  devices  and  developing 
particle-free  propellants  and  explosive  bolts.  Prelaunch  planning  of  future  space  operations  could 
include  the  use  of  programs  such  as  Collision  Avoidance  on  Launch  to  reduce  the  risk  of  collisions 
during  the  mission  [63] .  The  variability  of  launch  times  and  dates  are  limited  by  economic  feasibility 
and  mission  requirements.  Development  of  new  propellants  requires  time  and  resources.  However, 
securing  parts  to  the  spacecraft  is  easier  to  implement  and  thus  has  been  applied  in  some  cases  [49]. 

To  minimize  chances  for  explosions  of  discarded  rocket  bodies,  an  easy,  inexpensive,  and 
effective  technique  is  to  vent  residual  fuel  and  pressurants.  This  is  already  widely  practiced,  but 
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it  is  not  an  established  policy  by  all  space  agencies  or  companies.  Furthermore,  there  is  already  a 
considerable  number  of  propulsive  stages  left  in  orbit  by  previous  missions  with  residual  fuel,  and 
these  continue  to  present  an  explosion  hazard  as  illustrated  by  some  explosions  that  have  occurred 
over  a  decade  after  the  mission  was  launched  [31].  Also,  any  object  left  in  orbit  has  the  potential 
to  be  fragmented  by  collision. 

Inactive  payloads  form  a  reservoir  of  possible  sources  for  fragmentation  debris.  Disposal 
or  removal  of  satellites  at  EOL  is  more  important  for  GEO  where  orbit  lifetimes  extend  beyond 
thousands  of  years.  There  are  several  options  for  reorbiting  from  GEO: 

•  Deorbit 

•  Escape  trajectory 

•  Stable  longitudes 

•  Stable  plane 

•  Lower  altitude 

•  Raise  altitude 

Deorbiting  GEO  spacecraft  by  current  propulsion  technology  poses  unfavorable  economic 
penalties.  In  fact,  for  circular  orbits  above  25,000  km,  Earth-escape  maneuvers  are  less  expensive 
than  deorbiting  maneuvers  [51].  Passive  devices  may  be  deployed  or  inflated  to  increase  the  area- 
to-mass  ratio  which  increases  the  rate  of  orbital  decay,  but  this  raises  technical  difficulties  and  also 
increases  probabilities  of  collision  from  the  enlarged  satellite  cross-sectional  area. 

One  permanent  solution  is  for  active  satellites  to  make  a  propulsive  maneuver  at  EOL  to  an 
Earth-escape  trajectory.  However,  this  is  a  very  costly  option.  An  alternative  to  propulsive  burns 
to  escape  Earth  orbit  is  the  use  of  solar  sails  to  raise  the  satellite’s  altitude.  There  are  technical 
difficulties  in  the  deployment  and  control  of  these  solar  sails,  however,  and  the  technology  has  not 
been  tested  in  space.  Electric  propulsion  also  has  been  studied  as  a  feasible  means  to  boost  GEO 
satellites  into  a  heliocentric  orbit  [57]. 
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In  the  GSO,  there  are  two  stable  points,  one  at  75.3°E,  which  is  over  the  Indian  Ocean  and 
south  of  Bombay,  India,  and  the  other  at  255.3°E,  which  is  over  the  Pacific  Ocean  and  south  of 
Denver,  Colorado.  Because  uncontrolled  GEO  satellites  will  migrate  toward  the  closest  stable  point 
and  continue  to  oscillate  about  it  with  periods  of  several  years,  there  is  some  longitudinal  bunching 
about  these  points.  A  relatively  inexpensive  disposal  option  is  to  move  satellites  at  EOL  to  either 
one  of  these  points.  However,  objects  must  be  placed  with  negligible  position  and  velocity  errors; 
small  errors  in  the  rendezvous  maneuver  can  result  in  large  oscillations  around  the  stable  points  [13] . 
The  situation  is  further  complicated  by  the  presence  of  active  satellites  in  the  near  vicinity  of  these 
points.  A  collision  between  dead  payloads  at  the  stable  points  could  scatter  fragments  throughout 
GEO.  Most  researchers  do  not  favor  this  option  due  to  the  drawbacks. 

For  GEO,  a  “stable”  plane  exists,  inclined  7.4°  to  the  equator  [28].  The  inclination  of 
satellites  launched  into  this  inclined  orbit  remains  within  about  1.2°  without  any  stationkeeping. 
Therefore,  relative  velocity  between  satellites  in  the  stable  plane  is  small  compared  to  the  average 
relative  velocity  between  random  GEO  satellites  with  inclinations  as  low  as  0°  and  as  high  as  15°. 
The  wobbling  of  the  stable  plane  can  be  further  decreased  by  placing  satellites  at  specific  ascending 
nodes  [33].  Since  plane  change  maneuvers  are  generally  costly,  one  suggestion  is  to  place  satellites 
into  the  stable  plane  at  the  beginning  of  their  operational  lifetime.  In  doing  so,  however,  the 
advantages  of  geostationary,  as  opposed  to  geosynchronous,  orbits  would  be  lost. 

Subsynchronous  disposal  orbits  are  another  option.  The  probability  of  collision  decreases 
as  a  function  of  range  away  from  GEO;  it  is  two  orders  of  magnitude  less  at  100  km  below  GEO 
altitude  [13].  The  major  objection  to  this  solution  is  that  disposal  orbits  below  GEO  might  com¬ 
plicate  future  launches  as  new  payloads  pass  through  these  subsynchronous  orbits  on  their  way  to 
GEO. 

A  cost-effective  alternative  that  has  been  practiced  by  several  GEO  users  is  to  remove  the 
satellites  at  EOL  to  supersynchronous  storage  orbits,  also  called  disposal  or  graveyard  orbits.  To 
place  an  average-sized  GEO  spacecraft  about  300  km  above  GEO  altitude  requires  only  5-10  kg  of 
propellant  which  translates  to  a  1-2  months  reduction  in  the  satellite’s  mission  lifetime  [13].  The 
probability  of  collision  is  then  reduced  by  two-  or  three-orders  of  magnitude  [42].  This  solution 
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too  has  difficulties.  As  the  population  increases  in  the  disposal  orbits,  the  probability  of  collisions 
between  inactive  payloads  increases.  Fragments  from  collisions  in  storage  orbits  may  be  thrown  into 
orbits  that  intersect  GEO.  However,  this  is  the  only  economically  feasible  option  with  acceptable 
drawbacks  currently  available.  Despite  this,  some  feel  that  disposal  orbits  are  only  an  intermediate 
solution  [46,  56]. 

These  are  self-disposal  options,  and  they  are  not  without  risks.  There  is  some  danger  of 
malfunction  in  all  propulsive  systems  and  more  so  for  systems  that  already  have  been  subjected 
to  the  harsh,  natural  environment  of  space  for  several  years.  Further,  propellant  gauging  is  often 
inaccurate,  leaving  some  satellites  at  EOL  with  insufficient  fuel  for  the  disposal  maneuvers.  It  has 
been  noted  that  a  partially  successful  disposal  maneuver  may  create  a  worse  situation  than  none  at 
all  [63].  However,  to  remove  inactive  GEO  satellites  that  were  not  reorbited,  active  retrieval  with 
an  orbital  maneuvering  vehicle  or  by  attaching  a  propulsion  system  or  a  passive  device  to  increase 
area-to-mass  ratio  is  necessary.  The  technology  required  for  active  retrieval  in  GEO  is  available, 
but  it  is  not  cost-effective  at  the  present  time. 

1.5  Summary 

The  orbital  debris  situation  in  GEO  is  of  particular  concern  to  some.  Many  feel  that  keeping 
collision  probabilities  low  by  not  augmenting  the  drifting  GEO  population  is  of  great  importance. 
Besides  the  sheer  number  of  GEO  objects,  their  distribution  is  also  a  factor.  Due  to  orbital 
perturbations,  the  effect  of  orbital  concentration,  or  bunching,  at  certain  longitudes  is  already 
significant  [14].  Satellite  retrieval  or  disposal  after  their  active  lifetime  is  presently  not  practicable 
in  GEO.  Therefore,  active  EOL  disposal  schemes  are  the  only  economically  feasible  option  at  the 
present,  but  none  are  consistently  practiced  by  all  launching  groups. 
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Chapter  2 


Research  Overview 


Since  there  are  no  natural  sinks  for  debris  in  GEO,  the  accumulation  of  inactive  satellites 
and  other  operational  debris  in  this  valuable  regime  has  raised  particular  concern  among  the  many 
and  diverse  GEO  satellite  operators.  The  practice  of  raising  the  altitude  of  satellites,  placing  them 
out  of  GEO  before  stationkeeping  fuel  is  expected  to  be  depleted,  is  mainly  for  the  purpose  of  main¬ 
taining  continuous  operation,  but  it  is  also  a  cost-effective  mitigation  measure  that  many  operators 
are  voluntarily  implementing.  One  immediate  drawback  is  that  the  usage  of  supersynchronous 
orbits  (SSO)  as  a  storage  area  for  debris  may  preclude  future  plans  to  utilize  it  in  a  more  produc¬ 
tive  manner.  However,  future  technology  may  allow  for  different  and  more  permanent  mitigation 
measures  for  GEO  and  perhaps  even  for  active  removal  of  debris  from  GEO  and  SSO,  if  the  need 
arises. 


2.1  Identification  of  Problem 

In  considering  the  extensive  use  of  SSO  for  satellite  disposal,  the  first  concern  is  to  determine 
the  minimum  safe  distance  above  GEO  such  that  objects  in  the  disposal  orbits  will  not  interfere  with 
the  GEO  population  in  the  future.  This  involves  both  defining  the  useful  GEO  area  and  studying 
the  perturbation  effects  on  objects  in  SSO.  Since  there  are  numerous  uncontrolled  objects  in  GEO 
as  well,  perturbation  analysis  for  these  objects  is  also  relevant.  Various  reorbiting  guidelines  have 
been  proposed  by  different  spacefaring  nations  and  agencies. 
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Previous  studies  have  focused  on  propagating  the  orbits  of  intact  objects  in  GEO  and  SSO. 
However,  in  the  aftermath  of  a  collision  in  SSO,  pieces  of  varying  sizes  and  shapes  can  be  found 
in  orbits  quite  different  from  the  parent  objects’  orbits.  Depending  on  the  altitude  of  the  parent 
objects  in  the  storage  orbits  before  the  breakup,  these  pieces  may  intersect  GEO  by  the  impact  of 
the  collision,  or  subsequent  perturbing  forces  over  time  may  cause  the  pieces  to  mingle  with  the 
GEO  population. 

This  research  examines  the  aftermath  of  collisions  in  SSO,  and  in  particular,  the  orbits  of 
the  fragmentation  debris.  A  typical  reorbiting  altitude  is  selected  to  determine  the  semimajor  axis 
of  the  parent  objects’  orbits,  assumed  circular.  Since  collisions  produce  a  wide  variety  of  debris 
sizes  and  shapes,  an  average  piece  is  chosen  such  that  studies  of  the  orbit  of  this  one  representative 
piece  will  provide  insight  into  the  relative  spatial  distributions  of  smaller  and  larger  fragmentation 
debris.  Searching  for  a  worst  case  scenario,  the  orbit  of  this  fragment  is  analyzed  to  determine  if  at 
some  point  after  the  collision  it  will  drift  into  some  part  of  the  GEO  region.  Reiterating  this  process 
with  varying  reorbiting  altitudes  will  determine  a  minimum  altitude  for  the  storage  orbits  such  that 
an  unacceptable  fraction  of  the  collision  fragments  will  not  interfere  with  the  GEO  population. 

2.2  Earlier  Studies 

The  effects  of  orbit  perturbations  in  the  GEO  region  have  been  well-studied  for  many  pur¬ 
poses.  A  main  concern  of  satellite  operators  is  determining  the  necessary  stationkeeping  maneuvers 
and  the  amount  of  fuel  required.  For  orbital  debris  analysis,  the  motion  of  uncontrolled  satellites  is 
important  in  estimating  collision  hazards  and  predicting  trends  based  on  the  rates  of  orbital  decay, 
key  components  of  environmental  models. 

By  the  early  1980s,  various  GEO  operators  had  established  policies  for  reorbiting  at 
EOL  [61].  These  reorbiting  guidelines,  however,  differ  among  the  satellite  operators,  and  in  some 
cases,  they  differ  significantly  [13,  60].  Other  policies  do  not  specify  the  reorbiting  altitude  but  in¬ 
stead  prescribe  conditions  upon  which  the  reorbiting  altitude  may  be  calculated  as  appropriate  for 
individual  satellites.  The  Russian  Space  Agency,  for  instance,  requires  only  that  reorbited  space¬ 
craft  be  placed  such  that  it  will  not  approach  GEO  any  closer  than  200  km  [1].  As  more  studies 
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are  conducted  and  published,  movement  toward  consensus  on  reorbiting  altitudes,  given  either  as 
a  fixed  number  for  all  cases  or  as  an  expression  with  variables  that  are  influenced  by  individual 
spacecraft  characteristics,  may  be  possible.  Many  now  support  300  km  as  the  minimum  reorbiting 
distance  [46,  13].  Three  agencies  have  published  detailed  studies  leading  to  their  guidelines  which 
are  summarized  in  this  section.  These  studies  mostly  have  focused  on  defining  the  useful  GEO 
area  and  then  finding  the  minimum  altitude  above  it  such  that  debris  in  the  disposal  orbit  will  not 
likely  interact  with  the  GEO  population. 

2.2.1  NASA  Study 

The  NASA  Safety  Standard  [45]  defines  the  GEO  region  as  300  km  above  and  below 
35,788  km  altitude.  The  region  above  36,088  km  altitude  is  designated  as  the  disposal  region 
for  GEO  satellites.  The  Safety  Standard  and  a  study  by  Loftus  [41]  conclude  that  the  perigee 
altitude  in  the  storage  orbit  should  be  at  least 

GEO  +  300  km  +  2, 000  km  kg/m2  x  A/m  (2.1) 

with  the  area-to-mass  ratio  in  m2/kg.  The  first  term  accounts  for  perturbations  at  GEO  (50  km) 
and  at  the  disposal  orbit  (50  km),  operational  excursions  (50  km),  imperfect  insertion  at  the 
disposal  orbit  (50  km),  and  a  safety  margin  (100  km).  The  second  term  is  the  varying  effects  of 
solar  radiation  pressure  (SRP)  which  depends  on  the  objects’  area-to-mass  ratio. 

2.2.2  NASDA  Study 

Currently,  NASDA  requires  reorbiting  by  150  km  beyond  GEO  but  targets  a  desired  reorbit¬ 
ing  distance  of  500  km  [35,  36].  However,  perturbation  studies  [34]  have  led  to  working  expressions 
for  reorbiting  distances  based  on  individual  satellite  characteristics: 

A  a  =  186  km  +  0.011  kg/m2  x  a  x  Cr  x  A/m  (2.2) 


where 

A  a  =  reorbiting  distance  (km) 
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Table  2.1:  GEO  Perturbation  Effects  From  NASDA  Study. 


Perturbation  Source 

Radial  Variation 

Earth  oblateness 

35  km 

Solar  third  body 

53  km 

SRP 

45  km 

Total 

133  km 

a  =  semimajor  axis  of  the  disposal  orbit  (km) 

Cr  =  radiation  pressure  coefficient  of  the  satellite 
A  =  effective  cross-sectional  area  (m2) 
m  =  satellite  mass  (kg) 

This  equation  is  explained  as  follows.  The  186  km  in  Equation  2.2  is  the  sum  of  the  radial  variation 
that  a  GEO  satellite  experiences  due  to  all  perturbations,  as  detailed  in  Table  2.1,  and  the  radial 
variation  for  a  reorbited  satellite  in  SSO  due  to  solar  gravitation  (53  km).  Also  in  SSO,  Earth 
oblateness  is  a  much  smaller  effect  because  the  resonance  between  the  rotation  of  the  Earth  and 
the  mean  motion  of  the  satellite  is  shallow;  thus,  the  NASDA  study  allots  0  km  in  radial  variation 
for  Earth  oblateness  in  SSO.  The  entire  second  term  accounts  for  the  effects  of  SRP  on  a  reorbited 
satellite.  The  values  in  Table  2.1  were  computed  using  the  following  conditions: 

e  (eccentricity)  =  0.001 

A/m  =  0.005  m2/kg 
Cr  -  1.5 

This  equation  has  since  been  revised  [35,  36].  The  new  expression  is 

A  a  =  200  km  +  0.022  kg/m2  x  a  x  Cr  x  A/m  (2.3) 

While  reasons  for  the  14  km  increase  in  the  first  term  were  not  explicitly  stated,  variations  in 
assumed  satellite  orbital  elements,  namely  semimajor  axis  and  eccentricity,  will  explain  it.  The 
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additional  factor  of  2  in  the  second  term  is  to  account  for  the  long  term  variation  in  the  eccentricity 
vector.  The  original  expression  only  accounted  for  the  yearly  variation  in  eccentricity  vector. 

2.2.3  ESA  Study 

An  ESA  study  [5]  suggests  radial  limits  of  50  km  above  and  below  geosynchronous  altitude 
to  define  the  boundaries  of  GEO.  This  definition  is  based  on  the  actual  radial  variations  of  controlled 
satellites  in  GEO.  The  initial  formula  given  for  calculating  reorbiting  distance  to  null  the  collision 
rate  over  20  years,  as  explained  further  in  Reference  [5],  is 

da  =  1271  kg/m2Ae  +  56  km  (2.4) 

in  which  da  is  the  minimum  reorbiting  distance  above  GEO  in  km  and  Ae  is  the  effective  area-to- 
mass  ratio  in  m2/kg.  The  formula  is  explained,  as  follows.  1271  Ae  accounts  for  radial  variations 
due  to  SRP  (±928  Ae  km)  and  the  coupling  effects  of  SRP  and  the  gravitational  forces  of  Earth 
and  the  Moon  on  the  motion  of  the  eccentricity  vector  (±343  Ae  km).  The  56  km  is  the  sum  of  the 
width  of  a  0.1°  GSO  ring  (37  km)  and  19  km  to  account  for  the  perturbing  effects  of  Earth  and  the 
Moon.  Recognizing  that  the  eccentricity  variation  might  exceed  the  aforementioned  bounds  after 
20  years,  a  final  formula  with  additional  buffers  is  given: 

da  =  1600  kg/m2Ae  +  65  km  (2.5) 

No  additional  explanation  is  given  to  account  for  this  revision. 

More  recently,  Flury  [22]  states  that  ESA  recommends  a  minimum  reorbiting  distance 
above  GEO  of  300  km.  This  number  is  primarily  based  on  perturbation  effects  in  SSO.  The 
small  perturbations,  Earth’s  gravity  ( J2 ,  etc.)  and  third  body  effects  (solar  and  lunar)  cause  a 
±30  km  variation  while  SRP  effects  vary  with  area-to-mass  ratio.  For  A/m  =  0.1  m2/kg,  a  very 
conservative  value,  the  total  radial  variation  due  to  all  perturbations  is  ±200  km.  Smaller  area- 
to-mass  ratios  result  in  smaller  variations.  Therefore,  300  km  is  deemed  to  be  enough  to  account 
for  all  perturbations  acting  on  most,  if  not  all,  GEO-type  spacecraft  and  includes  a  100  km  buffer 
zone.  The  buffer  covers  GEO  perturbations  for  uncontrolled  objects  and  the  operational  zone  for 
active  satellites. 
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Table  2.2:  Comparison  of  the  Reorbiting  Equations. 


Study 

NASA 

NASDA 

ESA 

Year 

1992 

1996 

1996 

GEO  Perturbations 


Geopotential 

35 

Third  Body 

53  (solar) 

SRP 

45 

GEO  Total 

50 

133 

100 

SSO  Perturbations 


Geopotential 

mm 

0 

30 

Third  Body 

mi 

53  (solar) 

SRP 

20 

11 

7 

SSO  Total 

70 

64 

37 

Other 

200 

14 

0 

Reorbiting  Distance 

320  km 

211  km 

137  km 

2.2.4  Remarks  and  Comparison 

To  compare  the  reorbiting  equations  with  ESA’s  flat  300  km  recommended  distance,  values 
for  some  quantities  are  assumed.  For  an  average  area-to-mass  ratio  of  0.01  m2/kg,  Loftus’s  formula, 
Equation  2.1,  gives  320  km  as  the  appropriate  distance  to  reorbit.  The  latest  NASDA  formula, 
Equation  2.3,  gives  about  211  km  for  the  same  area-to-mass  ratio,  Cr  =  1.2,  and  a  =  42, 377  km. 
Loftus’s  guideline  is  noticeably  more  conservative  but  also  closer  to  the  300  km  recommendation. 

Next,  compare  the  individual  components  of  the  reorbiting  equations  as  far  as  possible. 
In  particular,  each  perturbation  effect  is  separated  out  to  gain  further  insight  on  the  differences 
between  independent  studies.  Table  2.2  summarizes  the  studies  discussed  above  with  A/m  = 
0.01  m2/kg,  a  =  42,377  km,  and  Cr  =  1.2.  Distances  are  radial  and  in  km. 
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Table  2.2  is  missing  some  numbers  that  could  not  be  deduced  from  the  published  studies  so 
it  is  difficult  to  focus  on  specifics.  The  studies  differ  on  GEO  perturbations  by  83  km,  while  SSO 
totals  are  within  33  km.  It  is  possible  that  another  independent  study  will  help  in  understanding 
these  discrepancies.  Note  that  the  reorbiting  distances,  as  reflected  in  this  table,  are  based  on 
specific  area-to-mass  ratios.  Some  studies  prefer  to  use  an  extremely  conservative  value  and  then 
recommend  one  reorbiting  distance  for  all  satellites.  Others  leave  the  final  distance  up  to  the  satellite 
operators  to  choose  based  on  information  for  individual  satellites.  Finally,  note  that  objects  are 
assumed  to  remain  intact  in  these  studies. 

2.3  Approaches  and  Tools 

This  section  outlines  possible  methods  and  tools  to  study  the  effects  of  collisions  in  SSO. 
Although  this  research  employed  one  particular  method,  alternatives  are  described  to  provide  ideas 
for  refinement  of  this  work  in  the  future.  One  part  of  this  problem  involves  using  a  low-velocity 
collision  model  to  obtain  the  typical  fragment  characteristics  and  delta-velocities  imparted  to  the 
fragments.  Currently,  low-velocity  models  are  few,  and  those  in  existence  have  incorporated  little, 
if  any,  empirical  data,  which  is  even  more  scarce.  The  collision  model  provides  an  idea  of  the  initial 
orbits  of  the  fragments;  the  other  major  piece  of  the  puzzle  is  determining  the  evolution  of  these 
orbits  over  time.  One  approach  is  to  conduct  an  analytic  study  on  the  effects  of  perturbations  on 
the  breakup  pieces  to  gauge  the  extent,  if  any,  to  which  the  breakup  debris  will  inhabit  the  GEO 
region.  The  varying  parameter  is  the  altitude  at  which  the  breakup  that  produced  these  fragments 
occurred.  From  the  results,  a  minimum  reorbiting  distance  above  GEO  can  be  approximated.  The 
perturbation  studies  can  also  be  carried  out  numerically,  but  searching  for  the  worst-case  scenario 
involves  methodical  guessing  of  the  initial  orbit  orientation  since  the  breakup  model  would  at  most 
provide  only  information  on  semimajor  axis  and  eccentricity. 

2.3.1  Low- Velocity  Collision  Models 

An  important  tool  in  studying  this  problem  is  a  collision  model  for  the  SSO  regime.  In  par¬ 
ticular,  information  on  the  typical  characteristics  of  and  delta-velocities  imparted  to  the  fragments 
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produced  from  nonhypervelocity  collisions  is  desired.  The  direction  of  the  delta-velocities  is  often 
assigned  randomly;  consequently,  the  delta-velocity  can  be  applied  in  the  direction  that  produces  a 
new  orbit  with  minimal  radial  distance,  since  the  worst  possible  case  is  desired.  Unfortunately,  the 
vast  majority  of  research  in  impact  breakup  modeling  centers  on  hypervelocity  collisions  associated 
with  objects  in  LEO.  Most  ground  test  and  on-orbit  breakup  data  are  also  geared  for  the  hyper¬ 
velocity  regime.  The  first  decision  then  is  whether  to  develop  a  new  collision  model  for  use  in  this 
study  or  to  try  to  utilize  one  of  the  few  largely  unvalidated  models  available.  The  few  documented 
GEO  models  are  discussed  next. 

Long  term  evolution  programs  for  GEO  generally' include  nonhypervelocity  collision  models 
as  a  source  of  debris  generation  from  the  on-orbit  population.  Logically  then,  searching  for  appro¬ 
priate  collision  models  begins  with  GEO  environmental  models.  One  such  program,  developed  at 
the  University  of  Colorado,  is  ODESI,  which  includes  two  nonhypervelocity  collision  models  [43]. 
The  first  is  a  simplistic  spherical  model  that  is  statistically-based.  Assuming  that  users  have  good 
data  for  the  number  of  fragments  produced  and  the  delta-velocities  imparted,  the  model  breaks  the 
parent  object  up  in  equal-sized  pieces  and  randomly  generates  directions  for  the  delta-velocities. 
For  the  velocity  imparted  to  breakup  fragments,  this  is  usually  determined  using  on-orbit  breakup 
data.  That  is,  with  knowledge  of  the  orbits  of  the  breakup  pieces  and  also  of  the  parent  object’s 
orbit  before  impact,  delta-velocities  can  be  estimated.  However,  it  is  unlikely  that  a  large  number 
of  fragments  from  GEO  or  SSO  collisions  could  be  detected  or  tracked.  With  no  better  num¬ 
bers  readily  available,  the  study  assumed  delta-velocities  that  were  equal  to  the  relative  velocities 
between  the  colliding  objects.  The  second  model  is  a  finite  element  approach,  using  the  general 
purpose  MSC/NASTRAN  program  in  conjunction  with  MSC/DYTRAN  to  simulate  the  dynamic, 
nonlinear  behavior  of  spacecraft  material.  The  results  from  this  model  depend  heavily  on  the  grid 
size  employed;  fragments  produced  cannot  be  smaller  than  the  element  grid  size.  The  published 
study  includes  results  from  a  0.5  cm  grid  and  using  one  specific  GEO  satellite.  It  has  not  been 
determined  whether  a  smaller  grid  size  would  be  more  appropriate.  Moreover,  this  model  does  not 
provide  the  velocities  imparted  to  the  fragments. 
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Another  model  was  developed  at  the  NTT  Radio  Communication  Systems  Laboratory  in 
Yokosuka,  Japan.  Fragment  mass  distribution  may  be  chosen  to  be  either  uniform  or  following  a 
power  law,  similar  to  hypervelocity  cases.  Some  low  velocity  impact  tests  performed  suggests  that 
mass  distribution  is  well  modeled  by  the  power  law  distribution  [68].  The  velocity  distribution  can 
then  be  generated  numerically  using  a  model  based  on  the  principle  of  conservation  of  momentum 
and  energy  [69].  Experimental  results  on  the  velocity  distribution  are  pending. 

2.3.2  Analytic  Approach 

To  determine  a  suitable  storage  altitude  in  SSO,  two  questions  are  considered: 

•  What  is  the  maximum  altitude  traversed  by  typical  GEO  objects,  operational  as  well  as 
uncontrolled? 

•  What  is  the  minimum  altitude  reached  by  an  object  ejected  from  a  collision  in  the  storage 
orbit? 

The  first  concerns  the  definition  of  a  GEO  region,  the  protection  of  which  is  the  primary 
motivation.  The  answer  partly  depends  on  the  activities  of  operational  satellites  and  is  found  in 
some  the  earlier  studies  of  GEO  reorbiting  requirements. 

For  the  uncontrolled  objects,  study  begins  by  examining  the  equation  of  motion  for  a  two- 
body  system  with  perturbations: 

?-=-£r-+/;  (2.6) 

for  which 

r  =  second  derivative  of  position  vector  with  respect  to  time 
r  =  position  vector  of  the  satellite  with  origin  at  center  of  Earth 
H  =  gravitational  parameter  of  central  body 
«  Gm© 

fp  =  perturbing  acceleration  per  mass  unit 
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Perturbations  may  also  be  studied  from  a  n-body  approach,  but  anticipated  difficulties  are  many. 
Variation  of  parameters  (VOP)  is  the  usual  tool  selected  to  manage  Equation  2.6,  and  for  this, 
there  is  the  choice  of  classical  Lagrangian  VOP,  Gaussian  VOP,  and  Hamiltonian  VOP.  The  last 
requires  some  familiarity  with  Hamiltonian  mechanics.  Variation  of  coordinates,  also  known  as 
Encke’s  method,  is  another  option  [15]. 

The  perturbing  acceleration  is  the  sum  of  individual  accelerations  from  several  sources. 
Objects  in  GEO  and  SSO  are  heavily  influenced  by  the  following  perturbations: 

•  Earth  gravity  harmonics 

•  Third  body  gravitation,  solar  and  lunar 

•  Solar  radiation  pressure 

Earth  gravity  harmonics  are  the  terms  of  a  mathematical  expansion  representing  the  Earth’s 
deviations  from  a  perfect  sphere.  The  largest  terms  of  the  zonal  and  tesseral  harmonics  are  J2  and 
,/22 ,  respectively.  J2  arises  from  Earth  equatorial  oblateness  and  causes  secular  variations  in  the 
right  ascension  of  the  ascending  node  and  the  argument  of  perigee.  J22  is  associated  with  the 
ellipticity  of  the  Earth  equatorial  plane.  This  generates  the  long  term  (860-day)  oscillation  of 
geosynchronous  satellites  which  can  be  counteracted  by  periodic  stationkeeping  maneuvers.  While 
J2  effects  are  much  more  pronounced  than  lunar  and  solar  gravitational  perturbation  at  low  alti¬ 
tudes,  the  effects  of  each  are  about  even  in  GEO. 

Earth  satellites  are  affected  by  the  gravitational  attractions  of  other  celestial  bodies  in  the 
solar  system.  The  perturbation,  or  acceleration,  is  proportional  to  the  mass  of  the  perturbing  body 
and  inversely  proportional  to  the  distance  between  the  body  and  the  satellite.  Lunar  gravitational 
perturbations  are  significant  because  the  Moon,  though  small  in  mass  relative  to  nearby  planets, 
is  the  closest  body  to  the  Earth-satellite  system.  Conversely,  the  Sun’s  perturbing  effects  are 
important  due  to  its  great  mass.  The  net  result  is  a  precession  of  the  satellite’s  orbit  about  an  axis 
normal  to  the  perturbing  body’s  orbital  plane.  Lunar  and  solar  effects  are  significant  for  orbits  at 
high  altitudes  and  those  with  periods  greater  than  12  hours  [15]. 
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Solar  radiation  pressure  (SRP)  induces  the  yearly  variations  in  eccentricity  and  argument 
of  perigee  for  GEO  satellites  due  to  the  apparent  annual  motion  of  the  Sun  [64].  Eccentricity  may 
vary  from  0.001  to  0.004  in  six  months  [15].  The  effect  is  proportional  to  the  satellite’s  effective  area 
and  surface  reflectivity  and  is  inversely  proportional  to  satellite  mass.  Very  small  particles,  such  as 
micrometeoroids  and  solid  rocket  effluents,  with  larger  area-to-mass  ratios  are  strongly  influenced 
by  this  perturbation  force. 

Equations  of  motion  have  been  developed  for  the  gravitational  effects  of  Earth’s  asphericity 
and  nonuniformity  in  terms  of  the  classic  Keplerian  elements  in  the  works  of  Kaula  and  oth¬ 
ers  [37,  38].  However,  in  working  with  near-circular  and  low  inclination  orbits,  such  as  those  of 
most  objects  in  GEO  and  the  storage  orbits,  a  nonsingular  set  of  orbital  elements  may  be  more 
appropriate.  There  are  many  nonsingular  sets  available  [64],  but  the  equations  of  motion  and 
perturbing  accelerations  will  need  to  be  derived  in  the  new  coordinates.  Several  published  papers 
feature  development  of  the  desired  expressions  in  nonsingular  equinoctial  elements  [8,  10,  11]  and 
other  nonsingular  sets  [44]. 

Two  approaches  to  determine  the  maximum  periodic  effects  predicted  by  the  equations 
of  motion  are  described.  The  first  is  a  frequency-independent  method  in  which  the  amplitudes 
of  the  periodic  effects  from  each  perturbation  are  simply  summed  without  regard  to  the  various 
effects’  periodic  phasing.  Special  care  must  be  taken  to  understand  the  frequency  of  the  periodics 
when  combining  the  effects  to  properly  calculate  the  maximum  change  in  satellite  altitude.  This 
approach  is  the  more  conservative  of  the  two  proposed.  The  second  method  accounts  for  the  phasing 
differences  in  the  periodic  effects.  This  approach  attempts  to  maximize  the  radial  rate  of  change 
by  searching  over  the  orbit  orientation  elements.  A  symbolic  manipulator  software  package  can 
be  employed  to  help  with  the  mathematics.  Then,  based  upon  the  frequencies  and  phasing  of  the 
periodics,  the  maximum  change  in  radial  amplitude  can  be  deduced.  These  are  idealistic  ways  to 
deal  with  very  complicated  expressions.  In  practice,  methodologies  need  not  completely  disregard 
phasing,  but  at  the  opposite  end,  accounting  for  all  phasing  effects  may  be  impractical.  For  an 
initial  study,  something  between  most  desirable  and  sufficient  is  appropriate. 
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Once  the  maximum  altitude  of  the  GEO  region  and  the  minimum  altitude  of  objects  ejected 
from  SSO  collisions  are  determined,  the  storage  altitude  can  be  evaluated.  This  altitude  is  the  key- 
parameter  to  vary  in  an  iterative  manner  to  determine  the  minimum  safe  reorbiting  distance  above 
GEO. 

2.3.3  Numeric  Approach 

For  this  problem,  a  numeric  approach  can  be  undertaken  as  an  independent  study  or  in 
conjunction  with  an  analytic  part.  Discussion  of  the  latter  follows  first.  Validation  of  analytic 
results  can  be  achieved  by  conducting  a  numeric  analysis  using  a  proven  orbit  propagator.  While 
special  perturbation  methods  are  often  the  preferred  choice  for  high  accuracy  numerical  studies, 
they  may  not  be  the  best  choice  in  this  case  since  individual  periodic  effects  cannot  be  easily  isolated 
nor  broken  down  into  short-  and  long-periodic  components.  This  is  an  important  capability  so  that 
each  part  of  the  analytic  work  can  be  thoroughly  checked.  While  selecting  a  propagator,  this  should 
be  kept  in  mind. 

The  other  numerical  approach  can  be  performed  without  the  analytic  portion.  This  is  an 
important  option  should  analytic  approaches  fail  to  give  satisfactory  results.  In  this  strictly  numeric 
scenario,  a  genetic  search  algorithm  can  be  employed  with  an  orbit  propagator  to  determine  the 
minimum  radial  distance  of  fragments  resulting  from  collisions  in  the  storage  orbit.  The  search 
algorithm  propagates  a  set  of  elements  from  an  element  set  family  prescribed  by  the  breakup 
model  and  storage  altitude  to  find  the  orbit  orientation  that  produces  the  lowest  possible  altitude, 
i.e.,  to  find  the  worst  case  scenario.  A  concern  here  is  the  length  of  time  to  propagate  these  orbits. 
Physics  ensure  that  debris  with  moderate  area-to-mass  ratios  will  reside  in  the  GEO  region  and  SSO 
for  a  great  many  years.  However,  technology  advances  in  the  next  50  years  or  so  may  enable  the 
practice  of  different  mitigation  measures.  It  therefore  seems  unnecessary  to  look  beyond  this  time 
frame.  So  again,  special  perturbation  methods  may  not  be  the  ideal  choice  here  due  to  concerns 
over  computation  time  and  efficiency.  Semianalytic  propagators  should  be  considered  for  their 
advantage  in  this  area.  It  is  anticipated  that  this  strictly  numeric  approach  requires  significant 
computing  power.  Genetic  search  algorithm  may  need  to  be  implemented  with  propagators  using 
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high  performance  parallel  virtual  machine  architecture  [66].  This  could  be  a  desirable  option  if  a 
strictly  numerical  approach  is  pursued. 


2.4  Chosen  Methodology 

This  research  employs  an  analytic  approach  complemented  with  numerical  results  for  vali¬ 
dation.  The  analysis  uses  a  variation  of  parameters  formulation  in  nonsingular  equinoctial  elements, 
the  specifics  of  which  are  summarized  in  Chapter  3.  SRP  is  treated  like  a  conservative  force  with 
the  assumption  that  objects  are  continuously  sunlit.  Some  frequency-independent  assumptions  are 
made.  The  analytical  study  will  provide  greater  physical  insight  on  the  perturbation  sources  and 
effects  than  a  completely  numeric  approach.  It  will  also  be  more  general  so  that  it  can  be  applied  to 
a  variety  of  scenarios.  For  the  complementary  numerical  part,  the  Cowell  numerical  integrator  and 
the  semianalytic  propagator,  Draper  Semianalytic  Satellite  Theory  (DSST),  both  contained  in  the 
Draper  R  &  D  version  of  the  Goddard  Trajectory  Determination  System,  were  chosen.  The  Cowell 
propagator  is  useful  for  checking  short  periodic  variations.  For  longer  integration  runs,  the  DSST 
averaged  orbit  generator  is  employed  because  it  allows  for  specific  model  tailoring  with  accuracy 
comparable  to  special  perturbation  methods  [7].  The  equations  of  motion  in  DSST  are  decomposed 
into  short-  and  long-periodic  contributions  with  the  amount  of  force  modeling  configurable  at  run 
time  [18,  23].  Thus,  isolation  of  the  periodic  effects  is  possible. 

It  is  hoped  that  a  clearer  understanding  of  the  effects  of  perturbations  on  objects  in  GEO 
and  SSO  will  be  obtained.  The  end  objective  is  to  provide  more  information  for  the  selection  of 
storage  orbits.  The  uniqueness  of  this  study  lies  in  the  inclusion  of  SSO  breakup  modeling  as  a 
deciding  factor.  Preliminary  studies,  using  hypervelocity  breakup  models,  suggest  that  breakups 
in  currently  proposed  SSO  might  well  be  a  hazard  to  GEO  [39]. 
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Chapter  3 


Perturbation  Models 


The  perturbation  models  utilized  in  this  research  are  presented  in  this  chapter.  In  con- 
sideration  of  space,  relevance,  and  usefulness,  only  the  final,  working  expressions  for  Lagrange’s 
Planetary  Equations  and  the  various  disturbing  functions  are  included  here;  for  detailed  deriva¬ 
tions,  the  reader  is  referred  to  the  listed  references.  For  this  background  material,  the  primary 
reference  is  Danielson  [19],  which  summarizes  the  mathematical  groundwork  for  the  Draper  Semi- 
analytic  Satellite  Theory  (DSST).  A  large  deviation  from  DSST  is  that  this  study  does  not  use  the 
single-  and  double-averaging  techniques  employed  by  DSST.  Since  averaging  can  be  a  tricky  and 
arduous  task  and  it  was  deemed  unnecessary,  this  study  did  not  employed  averaging  techniques. 

The  choice  of  orbital  elements  was  made  based  on  the  class  of  orbits  to  be  examined:  nearly 
circular  with  low  inclinations  in  the  GEO  region.  Hence,  a  nonsingular  set  is  preferred,  and  of  the 
various  sets  of  nonsingular  elements,  equinoctial  elements  were  selected.  A  very  tangible  advantage 
of  this  particular  orbital  element  set  is  the  fact  that  the  perturbation  theory  based  on  equinoctial 
elements  has  been  well-developed,  tested,  widely  published,  and  implemented  in  the  semianalytic 
propagator  based  on  DSST. 

3.1  Equinoctial  Elements 

This  section  provides  a  brief  introduction  to  the  equinoctial  elements,  with  definitions  and 
their  relationship  to  the  familiar  Keplerian  elements.  The  equinoctial  elements  are  a,  h,  k,  p,  q,  and 
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A.  Associated  with  this  orbital  set  is  an  equinoctial  reference  frame  (f,  g,  w) ,  which  is  reproduced  in 
Figure  3.1  from  Reference  [19]  for  the  direct  equinoctial  elements.  Element  a  represents  semimajor 
axis,  same  as  in  the  Keplerian  set.  Elements  k  and  h  are  the  f  and  g  components,  respectively,  of 
the  eccentricity  vector,  while  q  and  p  are  the  f  and  g  components,  respectively,  of  the  ascending 
node  vector  in  the  equinoctial  reference  frame.  The  eccentricity  vector  points  to  periapse  from 
the  center  of  mass  of  the  central  body  and  has  magnitude  equal  to  the  eccentricity  of  the  satellite 
orbit.  Likewise,  the  ascending  node  vector  points  to  the  ascending  node  from  the  central  body, 
with  magnitude  depending  on  inclination.  The  final  orbital  element  is  the  mean  longitude  A. 

The  conversion  from  Keplerian  to  equinoctial  elements  is  given  by 


a  =  a 

h  =  esin(o;  +  /fi) 
k  =  ecos(u>  +  IQ) 


(3.1) 

(3.2) 

(3.3) 

(3.4) 
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q 


(3.5) 

(3.6) 


=  |tan  cos  ft 

A  =  M  us  1 Q 

where  I  is  the  retrograde  factor,  defined  by 

+1  for  direct  equinoctial  elements 
—  1  for  retrograde  equinoctial  elements 
It  is  useful  to  have  also  the  reverse  conversion.  Begin  by  computing  an  auxiliary  angle,  (: 

h 


I={ 


sinC  = 


cos£  = 


Vf?  +  k2 
k 


Then  proceed  with 


a  =  a 


i 

sin  ft 
cos  ft 
u 
M 


Vtf  +  k 2 

7 r  ^ arctan  \Jp2  +  q2 


P 


\Jp2  +  q2 

q 

■s/p2  +  q2 

C-/ft 


(3.7) 


(3.8) 

(3.9) 


=  A-C 


(3.10) 

(3.11) 

(3.12) 

(3.13) 

(3.14) 

(3.15) 

(3.16) 


3.2  Equinoctial  Reference  Frame 

From  Figure  3.1  on  page  25,  it  is  easy  to  see  the  main  defining  points  of  the  equinoctial 
reference  frame: 

•  f  and  g  are  in  the  satellite  orbital  plane. 

•  w  is  perpendicular  to  the  orbital  plane  and  parallel  to  the  angular  momentum  vector. 

•  The  right  ascension  of  the  ascending  node  is  equal  to  the  angle  between  f  and  the  ascending 
node. 
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Quantitatively,  the  components  of  the  equinoctial  frame  basis  vectors  (f ,  g,  w)  in  the  Earth- 
centered,  nonrotating  coordinate  system  ( x ,  y,  z )  are  given  by 


f  = 


1 


1  +p2  +  q2 


g  = 


1  +  p2  +  q2 


w  = 


l  +  p2  +  q2 


1  ~p2  +  q2 
2 pq 
-2  Ip 
2  Ipq 

(1  +  p2-q2)I 

2  q 

2 P 

-2  q 

(1  -P2-  q2)I 


(3.17) 


(3.18) 


(3.19) 


3.3  Direction  Cosines 


The  expressions  for  the  disturbing  forces  contain  direction  cosines  as  a  matter  of  convenience 
and  tidiness.  These  are  defined  by 

a  =  zb  •  f  (3.20) 

P  =  ZB-g  (3-21) 

7  =  zb  •  w  (3.22) 

where  zb  is  the  unit  vector  from  the  center  of  mass  of  the  central  body  to  its  geographic  north  pole 
in  the  geopotential  disturbing  function.  For  third  body  perturbations,  zb  points  from  the  central 
body  to  the  third  body.  In  the  case  of  solar  radiation  pressure,  for  which  a  disturbing  function 
does  exist  given  the  assumption  of  the  satellite  being  always  sunlit,  zb  points  to  the  Sun  from  the 
central  body.  Additionally,  the  direction  cosines  are  related  by 


a2  +  /?2  +  72  =  1 


(3.23) 
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3.4  Lagrange’s  Planetary  Equations 


Variation  of  parameters  is  a  method  to  study  the  disturbed  motion  of  two  bodies.  As 
developed  by  Lagrange,  perturbations  modify  two-body  motion  by 


dr 

dt 

dv  u, 

- 1 — _r 

dt  r3 


(3.24) 

(3.25) 


where  7Z  is  the  disturbing  potential  function.  Lagrange’s  Planetary  Equations  are  derived  from  the 
above  equations.  Battin  [3,  pages  476-484]  presents  a  rigorous  derivation  of  the  Planetary  Equations 
using  Keplerian  elements.  This  process,  when  repeated  for  equinoctial  elements,  remains  essentially 
unchanged,  aside  from  the  use  of  Poisson  brackets  instead  of  Lagrangian  brackets  [19]. 

Expressed  in  equinoctial  elements,  Lagrange’s  Planetary  Equations  are 


a 

h 

k 

P 
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A  dX 

Bdn  k  hB  dn 
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2  a  dn 
~A~da  A(1  +  B) 


,«7 

dn\ 


+  ^  dk  )  AB  Iq'B'iP'y) 


(3.26) 

(3.27) 

(3.28) 

(3.29) 

(3.30) 

(3.31) 


where 


A  = 

yfjla 

(3.32) 

B  = 

Vi  -  h2  -  &2 

(3.33) 

C  = 

l  +  p2  +  q2 

(3.34) 

and  the  cross-derivative  operator  is  defined  by 


3.5  Geopotential  Disturbing  Function 


The  familiar  expression  for  the  geopotential  is 

N  min(n,M)  v  „ 

n(r,<f>,ip)  =  ^2  (  — )  Pnm  (sin  (Gim  cos  mi)  +  Snm  sin  m^)  (3.36) 

r  n=2  TO=0  '  r  ' 

where 


/i  = 

r  = 


AT  = 

M  = 


<t>  = 


Cnm  i  Snm  — 


Ip  = 


gravitational  constant  of  central  body 

radial  distance  from  central  body  to  satellite 

maximum  degree 

maximum  order  (M  <  N) 

mean  equatorial  radius  of  central  body 

associated  Legendre  function  of  order  m  and  degree  n 

geocentric  latitude 

gravitational  coefficients 

geographic  longitude 


The  conversion  of  this  expression  into  equinoctial  elements  involves  writing  the  equation  in  complex 
form  and  two  Fourier  series  expansions,  one  in  true  longitude  and  the  other  in  mean  longitude. 
Danielson  [19]  goes  through  the  basic  derivations  steps.  Further  details  that  may  be  helpful  are 
contained  in  various  conferences  papers  [8,  54].  The  final  product  is 


where 


oo  M  N  N 

«  =  Hi  E  EE  E 

j=— oo  m= 0  s=—N  n=max(2,m,|s|) 


m 

{GJms  +  iHJms)(Cnm  -  iSnm)  exp[i(j\  -  n*0)]} 


jmprmpm  jg—n—ljSpvw 


r  ns  * 


(3.37) 


Re{z}  =  real  part  of  z 

I  =  retrograde  factor,  defined  by  Eq.3.7 


ym 
v  ns 


i=lL 


(n-}-s)?(n— s)! 


2”  (n-m)!(2f4)!(2ga)! 


if  n  —  s  is  even 


if  n  —  s  is  odd 


(3.38) 
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For  the  interested  reader,  Reference  [19]  contains  alternative  methods,  including  recurrence  rela¬ 
tions,  to  calculate  V™,  GJms ,  and  H3ms. 

Two  parts  of  the  disturbing  function  merit  further  discussion  here.  The  first  pertains  to 
the  calculation  of  the  kernels  KJn~l's  of  the  Hansen  coefficients  xjn~l's .  The  kernels,  functions 
of  eccentricity  or  h  and  fc,  are  defined  by 

K~n~1,s(e)  =  e-*s-jIX"n_1’s(e)  (3.45) 

They  may  be  obtained  in  three  ways  that  do  not  necessarily  involve  the  direct  computation  of 
the  Hansen  coefficients.  The  first  method  is  for  kernels  with  j  —  0  and  with  the  first  superscript 
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negative,  i.e.,  n  >  —1: 

0 


jv— n— 1  ,s  _ 

J\Q  — 


,1+2  s 


where 


if  n  =  5  >  0 
ifn  =  s  +  l>l 

[  [<2n  -  *)  V -  (»  -  *)XT ”+‘'*]  if  -  >  *  +  2  >  2 

1 


(3.46) 


X  = 


VI  -  w  -  ¥ 


Appropriate  initializations  are 


K°  =  1 
K1  =  -i 


(3.47) 

(3.48) 

(3.49) 


The  second  method  of  calculating  the  kernels  is  by  infinite  series  representation.  Note 
that  more  than  one  series  have  been  investigated  for  use  in  generating  the  kernels  via  the  Hansen 
coefficients  [53].  The  series  of  choice  is 

K~n~l,s  =  (i-h2-  k2)-n+?  f;  Y^~X(h2 + k2r  (s-so) 

where 


K—0 


=  modified  Newcomb  operators 
r)  =  max(j  —  s ,  0) 
l  =  max(s  —  j,  0) 

The  modified  Newcomb  operators  recommended  for  use  are  calculated  by 

4(p+a)Ytf  =  2(2£  -  C)  Y^  +  (£  -  0  Y<3£  -  2(2£  +  Q  Y^Z{ 
-  (£ +  0  Y%Zl  +  2(2p  +  2<r  +  2  +  3C)  FpC4,.-i 

which  are  initialized  by 

ro,o  —  x 

The  third  method  is  for  use  in  calculating  general  kernels: 


(3.51) 

(3.52) 


(3.53) 


(3.54) 


— n— l,s 


(3  —  tc)(1  —  n  +  s)(l  —  n  -  s) 


{(3  —  »)(1  —  n)(3  -  2n)K~n,s— 


(2  —  n)  1(3  -  n)(l  -  n)  +  ^  K~n+1’s+j2(  1  -  n)K~n+3's 


(3.55) 
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This  recurrence  relation  requires  values  of  four  kernels,  K~n,s,  K~n+l’s,  XJn+2's  and,  K~n+3’s,  for 
initialization.  These  starter  kernels  can  be  obtained  by  the  second  method  described  previously 
in  this  section.  All  kernels  for  this  study  were  generated  with  the  first  two  methods  since  the 
maximum  degree  examined  is  N  =  4.  The  third  method  is  useful  for  higher  degree  models. 

The  second  portion  of  the  disturbing  function  that  is  examined  more  closely  here  is  the 
calculation  of  the  Jacobi  polynomials.  Special  Jacobi  polynomials  with  m  =  0  are  calculated  by 

KUl)  =  2s  (^rt1  -  72)~§Pns(7)  (3.56) 

where  Pns( 7)  are  associated  Legendre  functions.  In  general,  the  Jacobi  polynomials  are  calculated 
by 


21(1+  v  +  w)(2l  +  v  +  w  -  2)Piw(j)  = 

(21  + v  +  w-  1)[(2 1  +  v  +  w)(2l  +  v  +  w  -  2)7  +  v2  -  tn2]Pr-i(7)  (3-57) 

-2(l  +  v-l)(l  +  w-l)(2l+v  +  itOP^T) 


and  are  initialized  by 


p<r  =  1 


pvw  _  Q 


(3.58) 

(3.59) 


3.6  Third  Body  Disturbing  Function 

Perturbations  caused  by  a  third  body,  point-mass  gravitational  field  are  modeled  by  the 
following  disturbing  function  [3,  19]: 

P3  rcosyA 


n(r,^,t)  =  \ |( 


|Rs-r| 


Rz 


(3.60) 


where 


/i3  =  gravitational  constant  of  third  body 

R3(t)  =  vector  from  central  body  to  third  body 

Rz  =  magnitude  of  R3 
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Third  Body 


Central  Body 


Satellite 


Figure  3.2:  Geometry  of  Third  Body  Perturbations. 


r  =  vector  from  central  body  to  satellite 

r  =  magnitude  of  r 

<p  =  angle  between  r  and  R3 


Figure  3.2  illustrates  the  geometry  of  this  situation. 

The  expansion  of  the  third  body  disturbing  functions  closely  follows  that  of  the  geopo¬ 
tential  [16,  19].  One  exception  is  the  introduction  of  direction  cosines  (a,  (3, 7)  which  reflect  the 
position  of  the  third  body  with  respect  to  the  satellite.  The  third  body  disturbing  function  in 
equinoctial  elements  is 

{  co  Af  Af  /  n  \  n  1 

£  (2  -  M  -5- )  VnsY™Qns[Cs(a,  13)  -  iSs(a,  ft)]  exp(ijA) 

i?3jil'oos=0n=maX(2,S)  Vi?3'  J 

(3.61) 

where 


N  = 

$0s  = 

Vns  = 

yns  _ 

*  n  ' 


maximum  power  of  parallax  factor 


Kronecker  delta 


("~g)! 


2”  mw 


if  n  —  s  is  even 
if  n  —  s  is  odd 


[fc  +  ih  sgn(s  -  j)],s  Kfs 


(3.62) 

(3.63) 
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sgn(a;)  =  sign  function,  defined  by  Eq.3.44 


«r 

Qns  (7) 

Cs(a,p),Ss(a,  0) 


=  kernels  of  Hansen  coefficients 

_  dsPn(j) 
djs 

=  functions  of  the  direction  cosines  a  and  P 


(3.64) 


Danielson  [19]  gives  recurrence  relations  for  the  Vns  and  Qns  functions  for  computational 
efficiency.  The  Yj ls  functions  deserve  a  few  extra  comments.  Once  again,  kernels  of  Hansen  coef¬ 
ficients  must  be  calculated,  and  this  is  accomplished  as  described  in  the  geopotential  section  with 


one  exception.  Kernels  with  j  =  0,  s  =£  0,  and  the  first  superscript  nonnegative  are  calculated  as 
follows: 


2s— 1  2,5—1 


rS— 1,5 


if  n  =  s  —  1>1 
if  n  =  s  >  1 


Ko  ={ 

.  2Z+1K0~1’S  ~  (ntn+iy)g0~2,S  if«>*+l>2 
in  which  x  is  defined  again  by  Equation  3.47.  Proper  initializations  are 


(3.65) 


K%’°  =  1  (3.66) 

Kq’1  =  -1  (3.67) 

The  kernels  with  j  =  0  and  s  =  0  are  calculated  along  with  general  kernels,  by  using  the 
infinite  series  representation: 

OO 

Kf  =  (1  -  h2  -  k2r+t  •£,  + *7”  (*-<») 

K=0 

where  Y^v>K+t  are  modified  Newcomb  operators  which  are  calculated,  as  before,  using  Equa¬ 
tion  3.53. 

The  Cs(a,/3)  and  Ss(a,p)  functions  are  defined  by  [19,  page  32] 


Cs(a ,  0 )  +  iSs(a ,  0)  =  (a  +  ip)s  (3.69) 

and  are  calculated  by  the  following  recursion  formulae 

Cs+1(a,P)  =  aCs(a,p)~pSs{a,p) 

Ss+i(a,f3)  =  pC3{a,  P)  +  aSs(a,  p) 
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(3.70) 

(3.71) 


with  initializations 


Co  =  1  (3.72) 

So  =  0  (3.73) 


3.7  Solar  Radiation  Pressure  Disturbing  Function 

If  the  satellite  is  assumed  to  be  always  sunlit,  the  perturbation  from  solar  radiation  pressure 
(SRP)  is  a  conservative  force.  Hence,  it  may  be  represented  by  a  disturbing  potential.  For  the  class 
of  satellites  under  consideration,  this  shadowless  model  seems  reasonable,  given  the  situation  that 
GEO  satellites  enter  Earth’s  shadow  only  around  equinox  and  then  for  only  a  little  over  an  hour 
during  each  pass.  The  disturbing  function  for  SRP  effects  is  given  in  References  [6,  19]: 

n  =  -  r 
|R©  -  rl 

for  which 
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radiation  pressure  coefficient  of  satellite 

cross-sectional  area  of  satellite 

mass  of  satellite 

mean  solar  flux  at  1  AU 

speed  of  light 

1  AU 

vector  from  satellite  to  Sun 


(3.75) 


This  disturbing  function,  when  expanded  in  equinoctial  elements,  is  identical  to  the  third 
body  disturbing  function,  with  two  exceptions: 

•  fj.3  is  replaced  by  — T  and 

•  n  in  the  summation  starts  at  1  instead  of  2. 
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Thus,  the  SRP  disturbing  function  in  equinoctial  elements  is 


TZ  =  Re 


{- 


T_ 
R © 


oo  AT  M 

E  E  E 

jzz—oo  s=0  n=max(l 


(2  -  50s)  Vnsy/sQns[Cs(a,  /?)  -  iSs(a, /3)]  exp(ijA)  > 

(3.76) 


All  functions  are  defined  as  in  the  previous  section. 
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Chapter  4 


Perturbation  Analysis  and  Results 


The  perturbation  models  detailed  in  Chapter  3  were  implemented  to  study  the  perturbed 
motion  of  uncontrolled  objects  in  geosynchronous  (GEO)  and  supersynchronous  (SSO)  regions. 
This  chapter  describes  the  methodology  employed  in  the  analytical  part  of  the  study,  validation 
of  the  models  by  numerical  integration,  and  analysis  of  each  perturbation  source  to  determine  its 
effects.  Application  to  the  reorbiting  problem  follows  in  Chapter  5. 


4.1  General  Methodology 

The  methodology  is  described  in  two  parts.  The  first  concerns  the  process  of  transforming 
the  expressions,  Lagrange’s  Planetary  Equations,  into  a  form  suitable  for  analysis.  This  involves 
obtaining  the  disturbing  potential  for  each  perturbation  source  in  workable  form,  taking  the  partial 
derivatives  of  the  potential  with  respect  to  the  orbital  elements,  and  substituting  these  expressions 
into  the  Planetary  Equations.  Some  assumptions  are  introduced  to  reduce  the  expressions  to 
manageable  sizes.  In  the  second  part,  the  variation  equations  are  examined  term  by  term  to 
determine  individual  effects  which  are  then  summed  to  calculate  the  cumulative  variations. 

4.1.1  Formulation  of  Lagrange’s  Planetary  Equations 

Evaluation  of  Lagrange’s  Planetary  Equations  3.26-3.31  can  be  a  cumbersome  job  if  con¬ 
ducted  manually.  In  this  study,  a  well-tested  and  widely  available  symbolic  manipulator,  Mathe- 


37 


matica,  is  utilized  to  handle  the  expressions  and  perform  somewhat  tedious  tasks,  such  as  expanding 
summations  and  taking  partial  derivatives  of  long  and  complicated  expressions.  The  great  advan¬ 
tage  of  this  route  is  that  common  errors,  such  as  dropped  terms  and  wrong  signs,  are  confidently 
avoided.  The  drawback  is  that  one  does  not  always  know  all  of  the  ramifications  associated  with 
built-in  functions  and  commands;  a  seemingly  straightforward  function  may  produce  quite  unex¬ 
pected  results  in  some  cases  but  not  others.  Blind  application  will  likely  produce  surprising  and 
erroneous  results;  therefore,  intermediate  expressions  should  be  examined  carefully. 

For  each  perturbing  force,  with  solar  third  body  and  lunar  third  body  treated  separately, 
the  associated  disturbing  potential  is  modeled  in  Mathematica.  Of  particular  concern  at  this  junc¬ 
ture  is  the  truncations  of  the  infinite  series  represented  by  summations.  In  general,  these  limits 
are  set  as  appropriate  for  the  orbit  class  under  consideration:  orbits  with  GEO-type  altitudes  and 
low  eccentricities.  Further  discussion  of  the  truncation  decisions  for  each  perturbing  force  follows 
in  later  sections.  The  orbits  of  GEO  and  SSO  objects  generally  have  inclinations  lower  than  15° 
which  dictates  that  direct  equinoctial  elements  are  proper;  thus,  the  retrograde  factor  I  is  equal 
to  1.  The  disturbing  potentials  are  composed  of  many  custom  functions  as  well  as  well-known 
ones,  such  as  Legendre  and  Jacobi  polynomials,  which  are  available  from  the  Mathematica  inter¬ 
nal  library.  The  custom  functions,  defined  and  verified  for  this  study,  are  fairly  straightforward, 
requiring  only  familiarity  with  Mathematica  and  general  programming  techniques.  After  the  sum¬ 
mations  are  expanded,  the  real  part  of  these  complex  expressions  form  the  disturbing  potential. 
Mathematica  performed  well,  except  for  one  small  matter.  The  factor  Vl  —  h~  —  k'2  appears  often 
in  the  disturbing  functions.  Not  knowing  that  h?  +  A:2  <  1  for  elliptical  orbits  so  that  the  factor  is 
always  real,  Mathematica  treats  it  as  possibly  complex,  thus  generating  some  extra  terms.  These 
terms  are  easily  identified,  and  discarding  them  properly  from  the  disturbing  potential  expressions 
is  a  simple  matter. 

Once  the  disturbing  potentials  are  available,  the  Planetary  Equations  can  be  formed.  This 
involves  taking  partial  derivatives  of  the  disturbing  functions  with  respect  to  orbital  elements  and 
direction  cosines,  which  Mathematica  handles  smoothly.  Since  the  study  is  primarily  concerned 
with  orbits  in  terms  of  their  radial  distance,  only  the  Planetary  Equations  addressing  changes  in 
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semimajor  axis,  Equation  3.26,  and  eccentricity,  i.e.,  equinoctial  elements  h  and  k,  Equations  3.27 
and  3.28,  are  of  interest.  Some  of  these  equations  for  a,  h ,  and  k  are  quite  large,  about  40, GOO- 
73, 000  lines  in  length  with  a  typical  visual  editor.  Therefore,  assumptions  are  made  to  reduce  the 
equations  to  a  workable  size.  The  first  is  a  zero  inclination  assumption  which  translates  to  p  =  q  =  0 
in  equinoctial  elements.  This  decision  is  based  on  the  fact  that  GEO  satellites  and  those  in  storage 
orbits  are  at  relatively  low  inclinations,  with  about  15°  being  the  upper  limit.  Discussion  of  the 
effects  of  this  assumption  is  contained  in  later  sections.  Mathematica  was  employed  to  implement 
this  assumption.  Next,  the  c,  h,  and  k  equations  are  reduced  to  first  order  in  eccentricity.  This 
works  out  to  be  roughly  equivalent  to  first  order  in  h  and  k,  also,  since  h  <  e,  k  <  e,  and  hk  <  e2. 
First  order  eccentricity  analysis  seems  appropriate  since  e  remains  relatively  small  for  uncontrolled 
objects  in  GEO  regime.  For  example,  a  GEO  satellite  with  initial  eccentricity  of  zero  and  subjected 
to  perturbations  (4x4  geopotential  field,  solar  and  lunar  gravitation,  and  solar  radiation  pressure) 
will  see  its  eccentricity  increase  to  about  0.0008  over  100  years,  according  to  one  study  [27].  Among 
the  reorbited  objects,  again  eccentricity  remains  small  [12].  The  orbits  of  collision  pieces  is  another 
matter.  However,  if  the  fragments’  orbits  have  eccentricities  larger  than,  say,  0.0035,  for  a  typical 
reorbiting  altitude,  the  fragments  are  likely  to  be  immediately  crossing  the  GEO  region.  Thus, 
perturbation  analysis  is  unnecessary  in  cases  of  large  eccentricity.  Removing  the  higher  order  terms 
from  the  a,  h,  and  k  equations  was  performed  manually  with  the  resulting  equations  fed  back  into 
Mathematica  to  be  expanded. 

4.1.2  Analysis  of  Lagrange’s  Planetary  Equations 

The  models  were  then  transferred  to  Microsoft  Excel  worksheets  for  term  by  term  analysis. 
At  this  point,  for  each  perturbation  source,  there  are  three  expressions,  one  each  for  a,  h,  and  k. 
In  order  to  estimate  the  amplitudes  of  the  variations  in  a,  h ,  and  k,  for  individual  terms  in  each 
expression,  analytical  integration  of  a  sort  is  performed.  This  is  generally  described  as  follows.  A 
typical  term  in  the  variation  equations  is  of  the  form 

a  =  o  cos  9  (4.1) 
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in  which  a  represents  one  term  in  either  the  equations  for  the  variations  in  semimajor  axis  or  in 
eccentricity.  Variables  with  no  secular  changes  are  integrated  as  constants.  Semimajor  axis,  eccen¬ 
tricity,  and  inclination  experience  no  secular  changes  in  the  GEO  regime  [12];  thus,  the  elements  a, 
h ,  and  k  are  treated  as  constants.  In  the  example  given  by  Equation  4.1,  a  is  a  quantity  subject 
only  to  periodic  variations.  Fast  variables,  9  in  the  example,  are  identified  and  approximated  by 
the  product  of  the  variables’  secular  rates  and  time: 

a  =  a  cos  9sect  (4.2) 

The  variables  describing  the  positions  of  the  satellite  and  the  perturbing  sources  within  their  orbits 
are  the  fast  variables.  Next,  integrate  over  time.  Equation  4.2  becomes 

a  a  •  n 

A  a  =  - —  sin  6  sect 

0$ec 

=  t— —  sin  9  (4.3) 

9  sec 

Initial  conditions,  90  in  this  example,  are  ignored  because  they  do  not  affect  the  amplitude  of  Ao 
and  because  this  particular  aspect  of  the  phasing  of  effects  is  disregarded. 

Two  points  concerning  this  analytical  integration  merits  further  discussion.  The  first  per¬ 
tains  to  the  amplitude  of  the  variations.  In  Equation  4.3,  the  amplitude  of  Aa  is  clearly  a/9sec. 
However,  although  the  integration  process  assumes  a  constant  “mean”  value  for  semimajor  axis 
and  eccentricity,  the  values  actually  employed  are  osculating  values.  Figure  4.1  illustrates  this 
difference.  Point  A  represents  the  worst  case,  in  which  the  osculating  value  is  Aa  larger  than  the 
mean  value.  From  this  point,  it  is  seen  that  the  actual  variation  is  twice  the  amplitude,  or  2Aa. 
Although  at  other  points,  the  variation  will  be  overestimated,  a  worst  case  bound  on  the  varia¬ 
tions  is  sought  here.  Therefore,  amplitudes  are  doubled  to  represent  the  largest  possible  variation. 
A  second  concern  with  this  analytical  integration  method  is  that  it  ignores  variations  caused  by 
changes  in  a,  h,  and  k.  For  semimajor  axis,  the  variations  in  a  are  very  small  in  proportion  to  the 
value  of  a  at  GEO  altitudes.  The  result  is  that  only  minor  errors  are  introduced  by  ignoring  the 
variations  in  a.  This  argument  does  not  work  for  eccentricity  or  h  and  k.  In  this  case,  both  eccen¬ 
tricity  and  the  variations  in  eccentricity  due  to  perturbations  are  small  and  roughly  of  the  same 
order  of  magnitude.  Fortunately,  the  most  significant  terms  in  the  a,  h,  and  k  equations,  i.e.,  the 
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Point  A 


Figure  4.1:  Variations  from  Osculating  Values. 


terms  causing  the  biggest  variations,  are  independent  of  h  and  k.  Thus,  these  terms  are  not  subject 
to  errors  caused  by  assuming  a  fixed  eccentricity.  The  eccentricity-dependent  terms  are  relatively 
small  and  do  not  seem  to  adversely  affect  the  calculated  maximum  variations  in  semimajor  axis 
and  eccentricity. 

Details  pertaining  to  each  perturbation  source  are  further  discussed  in  following  sections. 
Calculation  of  the  amplitudes  of  perturbation  effects  is  accomplished  using  Microsoft  Excel,  a 
spreadsheet  software.  Amplitudes  are  summed  by  unique  frequencies  and  by  all  frequencies.  The 
first  allows  for  the  identification  of  individual  effects  and  the  latter,  for  an  upperbound  on  cumulative 
effects,  since  the  phasing  of  the  individual  frequencies  are  disregarded.  As  explained  previously, 
the  amplitudes  are  doubled  to  approximate  actual  variations.  The  equinoctial  elements,  h  and 
fc,  are  related  back  to  Keplerian  eccentricity  for  physical  insight.  Examination  of  associated  h 
and  k  equations  shows  terms  with  similar  amplitudes  and  frequencies.  Furthermore,  variations 
in  eccentricity  are  contained  separately  and  fully  in  either  the  h  or  k  equations;  evaluation  of 
both  equations  is  been  redundant,  except  for  its  use  in  double-checking  details.  Thus,  for  each 
perturbation  model,  maximum  variations  in  semimajor  axis  and  eccentricity  are  generated  given 
initial  orbital  elements  and  the  object’s  physical  characteristics. 

The  results  are  validated  by  comparison  with  numerical  integration.  The  Draper  R  & 
D  version  of  the  Goddard  Trajectory  Determination  System  (GTDS)  contains  strictly  numerical 
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integrators  as  well  as  the  Draper  Semianalytic  Satellite  Theory  (DSST)  [24,  25,  48].  The  Cowell 
numerical  integrator  is  simple  to  set  up,  and  it  easily  handles  shorter  propagations.  Thus,  element 
history  plots  from  Cowell  integrations  provide  the  needed  information  on  short-periodic  variations. 
Long- periodic  variations  are  better  represented  by  the  DSST  average  orbit  generator  (DSST  AOG), 
in  which  short-periodic  terms  have  been  averaged  out.  The  DSST  integrations  in  this  study  are 
limited  to  20  years  because  of  the  limited  availability  of  accurate  third  body  ephemerides.  This 
is  deemed  sufficient  for  validation  of  the  analytical  models,  at  least  for  an  initial  study.  For  each 
perturbation  source,  the  significant  short-  and  long-periodic  variation  terms  in  semimajor  axis  and 
eccentricity  are  singled  out  and  compared  by  both  amplitude  and  frequency  using  the  same  test 
case.  Specifics  for  the  propagation  runs  are  given  in  Table  4.1  in  which  settings  are  the  same  for 
both  Cowell  and  DSST  AOG  unless  otherwise  noted.  Section  4.2.2  elaborates  on  the  choice  of 
semimajor  axis  and  eccentricity  for  the  test  orbit.  Inclination-dependent  effects  are  captured  by 
starting  the  integrations  at  15°  for  inclination;  for  other  variations,  inclination  is  initially  set  at 
zero.  Initial  values  for  node,  argument  of  perigee,  and  mean  anomaly  are  arbitrarily  chosen  to  be 
zero.  Note  that  the  400-second  step  size  for  the  Cowell  integrator  is  small  enough  to  allow  for  over 
200  steps/revolution.  Finally,  cumulative  perturbation  effects  from  the  analytical  and  numerical 
models  are  compared  in  three  cases  with  varying  semimajor  axes. 

4.2  Geopotential 

Particulars  of  the  methodology  concerning  only  the  geopotential  are  presented  first  in  this 
section.  This  includes  discussion  of  the  truncation  decisions  on  the  infinite  series  in  the  disturbing 
potential,  the  effects  of  the  assumptions  employed  in  reducing  the  variation  equations,  and  the  term 
by  term  analysis  by  integration.  Validation  of  the  geopotential  model  by  numerical  integration 
follows,  and  this  section  ends  by  summarizing  the  effects  on  semimajor  axis,  eccentricity,  and 
minimum  radial  distance  from  Earth’s  nonspherical  gravitational  forces. 
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Table  4.1:  Setting  Descriptions  for  Propagators  in  Test  Case. 


Integrator 

12th  order  summed  Cowell/ Adams  predict-partial  correct 

Step  Size 

400  s  for  Cowell 

43200  s  for  DSST  AOG 

Geopotential  Model 

4x4  JGM-2 

Third  Body  Model 

Solar/Lunar  based  on  JPL  DE  118/200  ephemerides 

SRP  Model 

spherical,  conical  shadow,  single  C r 

Semimajor  Axis 

42194  km 

Eccentricity 

0.007 

Inclination 

0,  15  degrees 

Ascending  Node 

0 

Arg  of  Perigee 

0 

Mean  Anomaly 

0 

Initial  Time 

midnight  on  180ct98  for  Cowell 

midnight  on  180ct88  for  DSST  AOG 

Final  Time 

Varies  with  period  of  perturbation 

C  R 

1.2 

A/M 

1  m2  /  100  kg 

Coordinate  System 

B1950.0 
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42205 


Time  (days) 

Figure  4.2:  Semimajor  Axis  History  for  4x4  and  8x8  Geopotential  at  GEO. 


4.2.1  Geopotential  Analysis 

A  4x4  Joint  Gravity  Model-2  is  employed,  and  for  the  GEO  regime,  this  is  deemed  suffi¬ 
cient  [28].  Figures  4.2  and  4.3  are  history  plots  of  semimajor  axis  and  eccentricity,  respectively,  for 
both  4x4  and  8x8  gravity  fields,  as  integrated  by  DSST  AOG  over  19  years.  Clearly,  the  variations 
in  semimajor  axis  and  eccentricity  are  very  similar  for  both  cases,  though  the  4x4  model  departs 
a  little  from  the  better  model  at  the  latter  part  of  the  integration,  as  seen  in  Figure  4.3.  This 
disparity  is  due  to  higher  order  terms  and  roughly  translates  to  altitude  differences  on  the  order  of 
hundreds  of  meters.  Since  perturbations  cause  overall  variations  on  the  order  of  tens  of  kilometers, 
the  error  introduced  by  this  truncation  of  the  gravity  model  is  acceptable.  Further  comparison  of 
the  two  cases  indicates  that  the  difference  in  mean  anomaly  does  drift  off.  However  the  position  of 
the  objects  at  any  specific  time,  as  given  by  mean  anomaly,  is  unimportant  for  the  purposes  of  this 
study.  Thus,  the  4x4  gravity  field  is  judged  to  be  appropriate. 
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Truncation  of  the  j  series  in  the  disturbing  potential  function,  Equation  3.37,  is  not  easily 
related  to  physical  attributes  of  the  orbits  on  first  examination.  It  is  clearly  related  to  eccentricity, 
but  the  summation  bounds  do  not  allow  one  to  pick  out  all  first  order  eccentricity  terms.  It  is  also 
related  to  frequency  by  generating  multiples  of  A.  Without  a  better  guess,  \j\max  is  initially  set  to 
2,  with  the  hope  that  all  significant  terms  would  be  captured.  However,  from  numerical  integration 
plots,  the  existence  of  a  sizable  resonance  term  for  j  =  —3, 3  is  concluded.  Accordingly,  the  bound 
of  the  j  series  is  extended  to  \j\max  —  3.  This,  together  with  N  =  M  =  4,  forms  the  limits  of  the 
basic  geopotential  model  studied. 

The  Planetary  Equations  were  formed  and  reduced  as  described  in  Section  4.1.1.  Incorpo¬ 
rated  into  the  variation  equations  are  the  direction  cosines  associated  with  geopotential  perturba¬ 
tions  which  are 


a  = 

P  = 

7  = 


2  V 

1  +  p2  +  q2 
2  q 

1  +P2  +  q2 
l-.p2-<?2 
1  +p2  +  q2 


(4.4) 


Reducing  the  expressions  to  first  order  in  eccentricity  did  not  eliminate  any  important  terms. 
However,  implementation  of  the  zero  inclination  assumption  did.  Comparison  with  numerical 
integration  shows  two  missing  inclination-dependent  terms  that  are  significant  for  low,  but  not 
zero,  inclinations.  They  are  the  J?  short-periodic  term  in  a  and  the  J3  long-periodic  effect  in  e 
which  were  removed  by  setting  p  =  q  =  0.  Since  GEO  and  SSO  objects  do  range  up  to  15°  in 
inclination,  inclusion  of  these  terms  is  deemed  important.  Identifying  a  missing  term  is  in  itself  a 
complicated  process.  First,  information  about  the  term’s  period  and  dependence  on  certain  orbital 
parameters  such  as  eccentricity  and  inclination  and  associated  gravitational  coefficient  is  gathered. 
Next,  the  n,  m,  s,  and  j  indices  in  the  analytic  perturbation  model  are  examined,  one  by  one,  to 
find  the  ones  that  will  generate  the  observed  characteristics  of  the  missing  term.  Once  identified, 
it  is  simple  enough  to  retrieve  them  with  Mathematica  and  add  them  to  the  basic  geopotential 
model. 
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After  the  a,  h,  and  k  expressions  are  transferred  from  Mathematica  to  Excel  spreadsheets, 
analytical  integration,  as  described  in  Section  4.1.2,  of  each  term  in  the  equations  is  carried  out. 
The  variations  in  semimajor  axis  and  eccentricity  are  then  given  by  expressions  for  each  term  which 
can  be  evaluated  for  a  given  orbit. 

The  values  employed  for  a,  h ,  and  k  are  partially  prescribed  by  the  semimajor  axis  and 
eccentricity  of  the  orbit  to  be  studied.  However,  h  and  k  also  depend  on  the  orbit’s  argument  of 
perigee  u  and  right  ascension  of  ascending  node  fi.  Both  are  left  as  variables  to  adjust  as  necessary 
to  find  the  largest  cumulative  variation  amplitudes.  However,  the  values  for  u  and  Q  have  little 
effect  on  the  overall  perturbing  effect  because  terms  containing  h  or  k  are  small. 

The  fast  variables  are  the  satellite’s  mean  longitude  and  Greenwich  sidereal  time  which  are 
approximated  by  the  product  of  a  constant  value  for  Asec  and  9,  respectively,  and  time.  The  secular 
rates  of  each  are  calculated  as  follows. 

The  secular  variation  of  A  is  approximated  by 


A  =  n  +  uj2 

,sec  “1“  &J2fsec  "h  ^©,sec  H“  ^ 0,sec  +  A SEP ,  sec 


in  which 


n 

&J2  ,$ec 
Qj2,sec 

^©,sec 
^0,sec 
A SRP ,  sec 


=  mean  motion  of  the  satellite 
=  secular  variation  of  u  due  solely  to  J<i 

=  secular  variation  of  fi  due  solely  to  J2 

=  secular  variation  of  A  due  solely  to  solar  third  body 

=  secular  variation  of  A  due  solely  to  lunar  third  body 

=  secular  variation  of  A  due  solely  to  SRP 


(4.5) 


The  secular  rates  of  ui  and  Q  are  based  on  only  the  effects  of  J2  because  it  is  by  far  the  dominant 
term  in  the  geopotential,  larger  by  several  orders  of  magnitude  than  the  next  largest  term.  The 
terms  A©,sec,  A0)Seo  and  A sRP,sec  do  not  include  the  mean  motion. 

Finally,  6  is  given  by 


X  _  361°  x  (tt/1800) 
86400 


rad/s 


(4.6) 
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The  dominant  term  in  Equation  4.5  is  of  course  the  mean  motion,  n.  However,  n  &  6  for  objects  in 
GEO  and  SSO.  Therefore,  when  j  =  m,  the  smaller  terms  in  Xsec  become  important  in  determining 
the  amplitudes  of  these  long-periodic,  resonant  effects. 

4.2.2  Geopotential  Model  Validation 

The  geopotential  model  is  validated  through  comparison  with  numerical  and  semianalytical 
integration  results.  Semimajor  axis  and  eccentricity  are  examined  separately,  and  for  each,  the 
significant  effects  can  be  classified  as  either  short-periodic  (SP)  or  long-periodic  (LP)  terms. 

One  case  is  employed  to  validate  each  of  the  analytical  perturbation  models  individually.  For 
this  study,  the  semimajor  axis,  eccentricity,  and  inclination  of  the  test  orbit  should  reflect  those  of 
objects  in  GEO  or  SSO  regions  or  of  fragments  generated  by  collisions  in  SSO.  In  fact,  semimajor 
axis  can  be  chosen  to  coincide  with  GEO,  but  near  GEO  altitudes,  the  analytical  geopotential 
model  notably  exaggerates  the  variations  due  to  the  resonant  terms.  This  is  because  the  analytical 
integration  holds  a  and  therefore,  Asec,  constant,  forcing  a  deeper  resonance  in  the  equations  than 
is  actually  experienced  by  uncontrolled  objects  in  osculating  orbits.  At  around  30  km  away  from 
GEO,  the  geopotential  model  recovers  accuracy.  Proposed  reorbiting  distances  are  no  less  than 
100  km  above  GEO;  thus,  this  limitation  should  not  affect  the  primary  objective  of  this  research, 
studying  the  effects  of  collisions  in  SSO.  The  problem  of  defining  the  GEO  region  will  be  handled 
in  Chapter  5. 

At  GEO+30  km,  the  resonances  are  shallow  enough  for  the  analytical  geopotential  model 
to  recover  accuracy.  Validation  is  performed  for  a  =  42194  km  and  e  =  0.007,  as  listed  in  Table  4.1, 
which  provides  details  of  the  numerical  and  semianalytical  integration  settings.  The  eccentricity  is 
larger  than  those  of  most  GEO  and  SSO  objects;  this  value  was  selected  in  anticipation  that  the 
orbits  of  collision  fragments  which  will  have  slightly  larger  eccentricities. 

Table  4.2  summarizes  the  significant  periodic  effects  for  this  test  orbit  as  given  by  the 
numerical  and  analytical  models.  The  amplitudes  listed  have  not  been  doubled  and  so  reflect  half 
variations,  as  explained  in  Section  4.1.2.  Most  terms  in  the  analytical  expressions  have  slightly 
larger  amplitudes  than  observed  amplitudes  from  plots  generated  by  numerical  and  semianalytical 
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Table  4.2:  Significant  Periodic  Effects  from  Geopotential  Perturbations. 


Numerical  Amplitude 

Analytical  Amplitude 

Period 

Semimajor  Axis 

SP 

115  m 

126  m 

12  hr  and  24  hr 

LP 

10  km 

11.9  km 

0.8  yr,  1.2  yr,  and  2.3  yr 

Eccentricity 

SP 

3.6  e-05 

3.7  e-05 

24  hr 

LP 

4.6  e-05 

3.8  e-05 

28  yr 

integration.  This  is  anticipated  from  employing  a  mostly  frequency-independent  method.  For 
this  study,  however,  establishing  a  bound  on  the  variations  in  semimajor  axis  and  eccentricity  is 
the  primary  objective  in  modeling  perturbation  effects.  While  comparison  of  individual  effects  is 
significant  to  ascertain  the  reliability  of  the  model,  irregularities  in  individual  terms  do  not  matter 
as  much  as  the  performance  of  the  entire  model.  It  is  also  important  to  note  that  though  the 
numerical  models  are  essentially  used  as  truth  models,  the  numerical  amplitudes  in  Table  4.2  are 
only  estimated  values,  which  contain  uncertainties  from  the  process  of  interpretation  from  plots. 
Discussion  of  individual  effects  follow. 

Figure  4.4  is  the  semimajor  axis  history  plot  showing  the  SP  variations.  The  test  case 
is  propagated  over  two  days  with  the  Cowell  integrator  and  initial  inclination  of  15°.  There  two 
periodic  effects,  one  with  a  12-hour  period  and  the  other,  a  24-hour  period,  captured  in  the  plot. 
From  the  analytical  model,  the  12-hour  variation,  which  is  inclination-dependent,  is  calculated  to 
be  almost  five  times  larger  than  the  24-hour  effect.  The  combined  amplitude  of  these  variations  is 
estimated  from  the  plot  to  be  115  km,  as  listed  in  Table  4.2,  but  the  individual  amplitudes  are  more 
difficult  to  assess.  Long-periodic  effects  in  semimajor  axis,  resulting  from  the  resonance  terms,  are 
seen  over  4  years  in  Figure  4.5.  The  analytical  geopotential  model  indicates  that  there  are  three 
resonance  terms  contributing  the  largest  variations  in  semimajor  axis.  In  the  plot,  only  two  can  be 
readily  identified:  the  effects  with  2.3-  and  1.2-year  periods.  The  third,  which  has  a  shorter  period 
of  0.8  years,  seems  to  be  masked.  According  to  the  analytical  model,  the  1.2-year  effect  is  by  far 
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Figure  4.4:  Geopotential  Short-Periodic  Variations  in  Semimajor  Axis. 
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Figure  4.5:  Geopotential  Long-Periodic  Variations  in  Semimajor  Axis. 
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7.01  E-03 


Time  (min) 

Figure  4.6:  Geopotential  Short-Periodic  Variations  in  Eccentricity. 

the  dominant  term  with  amplitude  about  10  times  greater  than  those  of  the  other  two  resonance 
terms. 

J2  is  the  force  behind  the  SP  variations  in  eccentricity,  shown  in  Figure  4.6.  The  plot 
clearly  shows  one  periodic  term  which  matches  the  analytical  term  in  amplitude  as  well  as  period. 
The  one  item  in  Table  4.2  that  may  be  of  concern  is  the  LP  variations  in  eccentricity  since  the 
analytical  model  appears  to  underestimating.  The  term  in  question  is  a  J3  effect  and  is  known  to  be 
dependent  on  inclination,  as  illustrated  by  the  following  plots.  The  one  in  Figure  4.7  corresponds 
to  a  zero  inclination  case,  in  which  the  J3  term  does  not  appear.  Multiple  effects  appear,  but  none 
are  significant  since  the  variations  are  on  the  order  of  10-7.  At  i  =  15°,  however,  which  is  the 
case  shown  in  Figure  4.8,  J3  has  a  substantial  effect.  The  variation  shown  in  this  plot  appears  to 
have  a  period  of  about  38  years,  though  only  roughly  half  of  the  period  is  captured  by  DSST  over 
20  years.  The  J3  term  in  the  analytical  model  has  a  period  of  28  years  which  corresponds  to  the 
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Figure  4.8:  Geopotential  Long-Periodic  Variations  in  Eccentricity  at  15°  Inclination. 

argument  of  perigee  rate  of  the  satellite.  This  difference  in  periods  may  be  due  to  the  fact  that 
the  propagator  included  only  geopotential  perturbations.  The  addition  of  third  body  and/or  SRP 
effects  may  influence  this  effect.  For  example,  inclination  history  plots  for  integration  with  only 
gravity  harmonics  in  addition  to  two-body  dynamics  show  that  inclination  will  exceed  16°  if  the 
initial  inclination  is  set  near  15°.  This  does  not  reflect  the  situation  in  reality  because  third  body 
perturbations  were  not  modeled.  Another  explanation  for  the  discrepancies  in  amplitude  as  well  as 
period  is  that  the  analytical  term  is  a  different  term  which  has  not  been  identified  and  isolated  in 
numerical  plots.  This  implies  that  the  effect  shown  in  Figure  4.8  is  missing  in  the  analytical  model. 
Future  work  should  include  further  investigation  into  the  LP  variations  in  eccentricity. 

4.2.3  Summary  of  Geopotential  Perturbation  Effects 

Figure  4.9  summarizes  the  variations  due  to  geopotential  perturbations  in  the  radial  direc¬ 
tion.  Data  for  this  plot  is  generated  by  the  analytical  model,  except  near  GEO  altitude  due  to 
the  limitations  of  the  analytical  model  in  that  vicinity.  The  near-GEO  data  is  supplied  by  DSST. 
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i=7° 
i=1  5° 


Figure  4.9:  Geopotential  Induced  Variations  in  Radial  Direction. 


Recall  that  total  variations  in  semimajor  axis  and  eccentricity  are  calculated  by  summing  the  indi¬ 
vidual  terms  of  different  frequencies  and  that  amplitudes,  such  as  those  in  Table  4.2,  are  doubled 
to  account  for  the  use  of  osculating  instead  of  mean  elements.  Let  A ageo  and  Aegeo  represent  the 
total  variations.  Then,  maximum  radial  excursion  is  defined  and  calculated  by 

Ar  =  ri  —  r2 

—  u(  1  e)  (u  Ac^ep^l  (e  -|-  Ae^goj) 

=  Aa5eo(l  -  e)  +  aAegeo  -  AageoAegeo  (4.7) 

The  altitude  above  GEO  is  given  by  a  —  42, 164  km  for  varying  semimajor  axes  while  eccentricity 
is  fixed  at  0.007,  the  value  used  in  the  validation  process.  In  Figure  4.9,  geopotential  perturbations 
are  shown  for  three  different  inclinations.  At  GEO,  the  70  km  in  radial  excursion  is  mostly  due  to 
the  resonant  effects,  which  decrease  sharply  as  altitude  increases  above  GEO.  At  GEO+200  km, 
effects  are  below  10  km  for  all  three  cases.  The  differences  in  radial  excursion  due  to  inclination 
are  overwhelmed  by  the  resonance  terms  until  around  GEO+60  km.  They  then  steadily  increase 
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Table  4.3:  Analytical  Expressions  for  Significant  Geopotential  Terms. 


Term 

Period 

Semimajor  Axis 

SP 

-12nxl2R^C2ocTh^{l-\-C)  2  [2pq  cos  2A  +  (p2  -  q2)  sin  2A] 

12  hr 

3/i1/2R|C2oa-5^2  (k  sin  A  —  h  cos  A) 

24  hr 

LP 

90V1'2 R%a-7'2  {S33  cos[3(A  -  6)]  -  C33sin[3(A  -  0)]} 

0.8  yr 

12 p^Rla-5'2  {S22  cos[2(A  -  0)]  -  C22  sin[2(A  -  0)]} 

1.2  yr 

-30^/2i4<r9/2  {S42  cos[2(A  -  0)]  -  C42sin[2(A  -  0)]} 

1.2  yr 

-3 y}l2R%a~7l2  {S31  cos(A  -  0)  -  C3i  sin(A  -  0)} 

2.3  yr 

Eccentricity 

SP 

-3/i1/2i?®C,2oa  7/2cosA 

24  hr 

LP 

1.5/jR$C3oa-6  sin  i  sin(2.5  sin2  i  —  1)  sin  ui 

28  yr 

until  approximately  GEO+200  km,  at  which  point  the  difference  between  the  0°  and  15°  cases  is 
about  2.5  km,  and  the  differences  appear  to  cease  growing. 

Table  4.3  lists  the  analytical  expressions  for  the  largest  terms  in  the  geopotential  model. 
The  periods  are  approximate,  and  they  correspond  to  the  validation  case  orbit  with  a  =  42, 194  km 
and  e  =  0.007.  All  terms  are  in  equinoctial  elements,  except  the  LP  eccentricity  term.  This 
is  presented  in  Keplerian  elements  because  searches  for  it  in  equinoctial  form  were  unsuccessful. 
Explanations  for  the  variables  in  the  analytical  expressions  are  given  in  Section  3.5  and  Equa¬ 
tions  3. 1-3.6  and  3.34.  The  variations  are  mostly  independent  of  eccentricity  with  the  exception  of 
the  small  J2  SP  semimajor  axis  term.  In  general,  larger  terms  have  no  dependence  on  eccentricity. 
The  dependence  on  inclination  of  the  J2  SP  and  J3  LP  effects  appears  in  Figure  4.9  and  produces 
about  3  km  difference  in  the  radial  direction.  Resonance  terms,  the  long  periodics  in  semimajor 
axis,  are  by  far  the  most  significant  factor,  and  they  are  not  functions  of  eccentricity  nor  inclination. 
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4.3  Solar  Third  Body 


Implementation  of  the  solar  third  body  perturbation  model  follows  a  similar  route  to  the 
geopotential  model.  The  potential  function  and  variation  equations  are  much  smaller  and  therefore, 
easier  to  manage.  A  major  difference  between  the  third  body  potential  equations  and  the  geopo¬ 
tential  analysis  is  the  introduction  of  new  direction  cosines  which  reflect  the  position  of  the  third 
body  with  respect  to  the  satellite.  Subtleties  in  those  direction  cosines  proved  to  be  nontrivial. 

4.3.1  Solar  Third  Body  Analysis 

The  third  body  disturbing  function,  Equation  3.61,  includes  two  infinite  series  which  require 
truncation  decisions.  For  solar  perturbations,  an  appropriate  truncation  of  the  parallax  factor 

/  \  jV 

f  is  M  =  3.  This  conclusion  is  based  on  the  truncation  algorithm  programmed  into  DSST 
which  is  summarized  in  References  [19]  and  [20].  The  j  series  is  truncated  by  =  2  based  on 

preliminary  examination  of  significant  frequencies  from  numerical  integration  element  history  plots 
and  also  on  the  order  of  eccentricity.  From  the  results,  it  appears  that  no  major  terms  were  left 
out  by  the  above  truncations. 

The  third  body  is  assumed  to  be  in  a  circular  orbit.  Therefore,  R3  is  constant  and  equal  to 
the  mean  distance  between  Earth  and  the  Sun  [59]: 

R3  =  1.0000010178  AU  x  149597870  km/AU  (4.8) 

This  is  a  reasonable  assumption  since  the  actual  eccentricity  of  Earth’s  orbit  about  the  Sun  is 
small,  about  0.0167  [64],  and  no  considerable  perturbing  effect  is  left  unmodeled  as  a  result.  The 
methodology  then  traces  the  same  pattern  as  before.  The  resulting  variation  equations  are  first 
order  in  eccentricity  with  zero  satellite  inclination  assumed. 

The  direction  cosines  employed  for  third  body  analysis  are 

a  =  cos  0,3  cos  O3  —  sin  CI3  cos  i3  sin  83 

/?  =  sin  £^3  cos  83  +  cos  cos  sin  63  (4.9) 

7  =  sin  13  sin  #3 
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where 


ft3  =  right  ascension  of  ascending  node  of  the  third  body 
03  =  argument  of  latitude  of  the  third  body 
=  argument  of  perigee  +  true  anomaly 
is  =  inclination  of  the  third  body 

These  expressions  already  reflect  the  assumption  of  zero  inclination  for  the  satellite;  this  explains 
why  the  satellite’s  inclination  and  ascending  node  do  not  appear  in  Equations  4.9.  For  the  Sun  [59], 

ft3  =  0,  by  definition  (4.10) 

i3  =  23°26’21.448”  (4.11) 

and  03  is  left  as  the  fast  variable  along  with  satellite  mean  longitude,  A.  During  the  model  validation 
process,  it  became  apparent  that  a  significant  effect  was  missing,  and  the  reason  for  its  absence 
can  be  attributed  to  the  above  direction  cosines,  Equations  4.9.  The  period  of  the  missing  term 
is  roughly  28  years,  about  half  of  the  period  of  the  satellite’s  node  rate.  As  perturbation  effects 

increase  the  inclination  of  the  satellite,  the  node  is  no  longer  undefined.  The  satellite  node  is  a 

secularly  but  slowly  varying  quantity.  Thus,  it  is  generally  masked  in  the  equations  by  the  faster 
variables,  except  in  a  few  terms.  The  missing  effect  is  recovered  by  isolating  the  suspected  term, 
in  this  case  n  =  $  =  2  and  j  —  0,  and  substituting  direction  cosines  with  the  node  of  the  satellite 
restored: 

a  =  cos(ft3  -  ft)  cos  03  -  sin(ft3  -  ft)  cos  £3  sin  03 

(3  =  sin(ft3  —  ft)  cos 03  +  cos(ft3  —  ft)  cos  i3  sin  03  (4.12) 

7  =  sin  i3  sin  03 

The  analysis  of  the  solar  third  body  model  again  involves  term  by  term  integration.  Unlike 
the  geopotential  expressions,  the  secularly  varying  variables  do  not  conveniently  appear  together 
in  the  argument  of  one  trigonometric  function  per  term.  Trigonometric  identities  are  employed 
to  expand  the  variation  expressions,  and  luckily,  Mathematica  has  the  capability  to  perform  this 
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tedious  process.  Values  are  needed  for  Xsec,  03, sec,  and  Qsec.  Calculation  of  Asec  is  covered  in 
Section  4.2.1.  For  solar  third  body,  0z<sec  is  given  by 

03  ,$ec  =  ^3,  sec  +  h  ,sec 

~  /3,sec 

*  3lJ?600  rad/s  (4'13) 

where  o?3 )Sec  is  ignored  because  it  is  extremely  small  [59].  The  secular  rate  of  /3  is  an  average 
rate  based  on  one  revolution  per  365.25  days.  The  satellite  secular  node  rate  is  the  sum  of  three 
components: 

^sec  =  ^J2,sec  "h  ^©,sec  f^0,sec  (4.14) 

in  which 

&J2,sec  =  secular  variation  of  due  solely  to  J2 

^©,sec  =  secular  variation  of  Q,  due  solely  to  solar  third  body 

Q,0}Sec  =  secular  variation  of  Q  due  solely  to  lunar  third  body 

Again,  since  J2  dominates  the  geopotential,  inclusion  of  the  effects  of  higher  zonals  is  unnecessary. 
The  secular  variation  due  to  SRP  is  also  ignored  in  Equation  4.14  because  it  is  negligible. 

4.3.2  Solar  Third  Body  Model  Validation 

Though  the  solar  model  is  sound  near  GEO  altitudes,  validation  is  carried  out  at 
GEO+30  km  to  preserve  uniformity.  Orbital  elements  and  numerical  integration  specifics  for  the 
test  case  are  summarized  in  Table  4.1.  The  significant  effects  of  the  solar  third  body  perturba¬ 
tions  are  listed  in  Table  4.4,  from  which  it  can  be  seen  that  the  agreement  between  analytical  and 

numerical  models  is  good.  Individual  effects  are  discussed  next. 

The  SP  variations  in  semimajor  axis  are  shown  in  Figure  4.10.  Short-periodic  effects  from 
solar  third  body  perturbations  are  generated  numerically  as  the  difference  between  a  solar  with  J2 
model  and  a  JVonly  model.  Hence,  the  values  for  the  semimajor  axis  are  sometimes  negative  and 
are  on  the  order  of  hundreds  of  meters.  The  12-hour  periodic  effect  is  clearly  seen  over  two  days.  Its 
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Table  4.4:  Significant  Periodic  Effects  from  Solar  Third  Body  Perturbations. 


Numerical  Amplitude  Analytical  Amplitude  Period 
Semimajor  Axis 

SP  470  m  473  m  12  hr 

LP  none 

Eccentricity 

SP  1.45  e-05  2  e-05  24  hr 

LP  3.3  e-05  3.3  e-05  6  mo 

1.5  e-04  1.5  e-04  28  yr 


Figure  4.10:  Solar  Third  Body  Short-Periodic  Variations  in  Semimajor  Axis 


Figure  4.11:  Solar  Third  Body  Short-Periodic  Variations  in  Eccentricity. 


amplitude  is  somewhat  larger  than  the  SP  effects  due  to  geopotential.  There  are  no  long-periodic 
terms  for  semimajor  axis. 

For  eccentricity,  the  SP  effects  are  captured  in  Figure  4.11.  Again,  some  values  for  eccen¬ 
tricity  are  negative  because  the  plot  shows  the  difference  between  two  models,  as  explained  in  the 
previous  paragraph.  The  primary  effect  shown  in  the  figure  has  a  period  of  about  24  hours.  The 
difference  between  analytical  and  numerical  amplitudes  in  Table  4.4  may  be  due  to  the  location 
of  the  third  body  during  numerical  integration  and  the  phasing  that  this  introduces  between  the 
different  components  of  the  terms  forming  the  total  amplitude.  Finally,  Figure  4.12  features  the 
LP  variations  in  eccentricity.  This  is  a  DSST  AOG  plot  of  a  solar  plus  J2  model.  The  plot  shows 
rapid  oscillations  on  top  of  a  larger  periodic  term.  The  periods  are  difficult  to  estimate,  but  es¬ 
timations  from  the  history  plot  do  correlate  fairly  well  with  the  information  from  the  analytical 
model.  The  oscillations  are  about  6  months  in  period  while  the  larger  variation  is  estimated  to  have 
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Eccentricity 


Time  (days) 


Figure  4.12:  Solar  Third  Body  Long-Periodic  Variations  in  Eccentricity 


- e=0.007 


e=0.002 


e=0.0007 


Figure  4.13:  Solar  Third  Body  Induced  Variations  in  Radial  Direction. 


a  28-year  period.  Future  work  should  examine  the  possibility  of  significant  inclination-dependent 
terms.  Some  eccentricity  plots  from  long  propagations  suggest  that  these  effects  are  likely  to  exist. 

4.3.3  Summary  of  Solar  Third  Body  Perturbation  Effects 

The  solar  third  body  effects  as  given  by  the  analytical  model  are  shown  in  Figure  4.13.  It  is 
evident  that  for  different  values  of  eccentricity,  solar  effects  vary  greatly.  The  effects  increase  linearly 
with  altitude  in  all  cases,  and  the  increase  is  greater  with  larger  eccentricities.  Radial  excursions 
are  under  20  km  in  the  range  of  eccentricity  and  semimajor  axis  that  this  study  anticipates. 

The  analytical  expressions  for  the  significant  terms  in  the  solar  third  body  model  are  listed 
in  Table  4.5.  Periods  of  the  terms  are  based  on  the  validation  case  and  are  shown  primarily  to 
relate  these  expressions  to  the  amplitudes  given  in  Table  4.4.  In  general,  for  large  eccentricities, 
beyond  the  range  of  Figure  4.13,  the  eccentricity-dependent  long-periodic  eccentricity  variations  are 
dominant.  The  short-periodic  effects  on  both  semimajor  axis  and  eccentricity  are  not  dependent 
on  eccentricity,  and  none  of  these  terms  are  affected  by  inclination. 
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Table  4.5:  Analytical  Expressions  for  Significant  Solar  Third  Body  Terms. 


Term 

Period 

Semimajor  Axis 

SP 

—1.4902  •  10-17a5/2sin  2A 

12  hr 

-1.7314  •  10-16a5/2sin[2(A  -  03)] 

12  hr 

LP 

none 

Eccentricity 

SP 

-3.5119  -  10_17a3/2  cos  A 

24  hr 

-1.3358  •  10-16a3/2  cos(A  -  2 03) 

24  hr 

-3.9659  •  10~18a3/2cos(A  +  203) 

24  hr 

LP 

1.9074  •  10“16a3/2  k  cos  203 

6  mo 

2.1602  •  10-16a3/2hsin  203 

6  mo 

2.9504*  10-20o3/2(hsin2f2  —  k  cos2£2) 

28  yr 
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4.4  Lunar  Third  Body 


Analysis  of  the  lunar  third  body  disturbing  function  closely  parallels  that  of  the  solar  third 
body  model.  However,  because  the  lunar  orbit  is  complicated  and  somewhat  eccentric,  greater  case 
must  be  taken  in  making  simplifying  assumptions  which  may  exclude  important  effects.  Addition¬ 
ally,  the  parallax  factor  is  larger  in  the  lunar  case,  thus  forcing  the  infinite  series  to  be  truncated  at 
a  larger  AT.  This  results  in  very  large  expressions  for  the  potential  function  and  variation  equations. 

4.4*1  Lunar  Third  Body  Analysis 

Like  the  solar  model,  there  are  two  infinite  series  in  the  lunar  disturbing  potential  for  which 
truncation  decisions  are  needed.  The  first  is  the  parallax  factor.  For  the  lunar  model,  Af  =  8. 
This  conclusion  is  based  on  DSST  built-in  truncation  algorithms  [19,  20].  The  j  series  is  again 
truncated  by  \j\max  =  2.  Following  the  method  employed  in  the  solar  case,  this  decision  rests  upon 
the  frequencies  appearing  in  the  element  history  plots  from  numerical  integration  and  the  order  of 
eccentricity  desired.  From  the  results,  it  seems  that  these  truncations  are  appropriate. 

The  Moon  is  assumed  to  be  in  a  circular  orbit.  Hence,  Rs  is  constant.  From  Reference  [55], 
the  mean  distance  between  Earth  and  the  Moon  is 

Rs  =  384400  km  (4.15) 

Unlike  in  the  solar  case,  this  assumption  is  generous,  given  that  the  actual  eccentricity  of  the  lunar 
orbit  is  about  0.0549  [64].  However,  handling  varying  values  of  Rs  analytically,  as  dictated  by 
an  eccentric  orbit,  is  a  difficult  matter.  Keeping  in  mind  that  significant  terms  can  be  retrieved 
later,  if  missing  from  the  model,  this  assumption  was  made.  It  unfortunately  led  to  the  omission  of 
the  most  significant  effect  in  eccentricity  variations.  This  particular  term  frustratingly  evaded  all 
searches  for  it  in  the  equinoctial  formulation  with  current  treatments  of  the  direction  cosines.  It 
was  instead  found  in  Kaula’s  third  body  model  [64]  which  is  in  Keplerian  orbital  elements  and  was 
added  to  this  lunar  model.  The  term  is  a  9.1-year  effect  due  to  the  lunar  argument  of  perigee  rate. 
The  zero  satellite  inclination  and  first  order  eccentricity  assumptions  do  not  add  further  difficulties. 
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The  direction  cosines  initially  employed  in  the  lunar  model  are  given  in  Equations  4.9.  The 
lunar  node  has  a  period  of  about  6798  days  [59],  much  smaller  than  the  periods  of  0 3  and  A,  and 
therefore,  is  hidden  in  most  terms.  The  inclination  of  the  Moon  with  respect  to  Earth’s  equatorial 
plane  varies  from  18.28°  to  28.28°  with  the  same  period  as  that  of  the  node  [2].  However,  the 
inclination  does  not  vary  secularly  and  thus  was  held  constant  through  the  analytical  integration. 
The  values  assigned  to  lunar  node  and  inclination,  45°  and  18.28°,  respectively,  maximizes  the 
perturbation  effects.  Once  more,  63  and  the  satellite  mean  longitude  A  are  the  fast  variables. 

Analytical  integration  of  the  lunar  third  body  model  follows  the  procedure  described  for 
the  solar  model.  The  value  for  lunar  63tSec  is  calculated  by  the  following  expression  [59]  which  is 
similar  to  Equation  4.13: 

@3, sec  = 


Again,  uz}Sec  is  ignored  because  its  period,  about  3232  days,  is  much  longer  than  the  true  anomaly, 
about  27.321661  days  [59]. 

4.4.2  Lunar  Third  Body  Model  Validation 

The  orbit  for  validation  purposes  is  set  once  more  at  semimajor  axis  equal  to  42194  km 
and  eccentricity  of  0.007.  Table  4.6  itemizes  the  largest  effects.  Numerical  integration  employed 
the  option  settings  detailed  in  Table  4.1.  Generally,  the  semimajor  axis  variations  correlate  much 
better  than  those  of  eccentricity.  Long  periodics  appear  only  in  eccentricity. 

Figures  4.14  and  4.15  show  the  short-periodic  variations  in  semimajor  axis  and  eccentricity, 
respectively.  These  plots,  like  those  featuring  SP  effects  in  the  solar  case,  depict  differences  between 
a  lunar  with  J2  model  and  a  J2-onlv  model.  Two  periodics  are  evident  in  the  semimajor  axis 
variations,  a  large  12-hour  effect  in  combination  with  a  smaller  24-hour  periodic.  The  eccentricity 
plot  shows  only  24-hour  terms.  The  discrepancy  in  the  amplitudes  of  the  SP  eccentricity  variations 
may  be  due  to  position  of  the  Moon,  as  it  represented  differently  by  numerical  and  analytical  models. 


^3, sec  “1“  y*3,sec 
/3,sec 

^  -rad/s 


2360591.51 


(4.16) 
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Table  4.6:  Significant  Periodic  Effects  from  Lunar  Third  Body  Perturbations. 


Numerical  Amplitude 

Analytical  Amplitude 

Period 

Semimajor  Axis 

SP 

1  km 

1.019  km 

12  hr  and  24  hr 

LP 

none 

Eccentricity 

SP  3.1  e-05 


4.5  e-05 


24  hr 


Difference  In  Eccentricity 


Figure  4.15:  Lunar  Third  Body  Short-Periodic  Variations  in  Eccentricity. 


Figure  4.16:  Lunar  Third  Body  Long-Periodic  Variations  in  Eccentricity. 


Numerical  integration  employs  accurate  third  body  ephermerides,  which  are  directly  dependent  on 
the  epochs  of  integrations,  while  the  analytical  model  fixes  the  third  body  in  the  orientation  that 
results  in  the  largest  radial  variations.  Hence,  the  analytical  model  is  expected  to  overestimate  and 
bound  the  maximum  possible  variations. 

The  lunar  analytical  model  is  somewhat  uncertain  in  the  LP  effects  on  eccentricity.  Fig¬ 
ure  4.16  is  a  DSST  AOG  plot  from  a  lunar  plus  J2  model.  The  9. 1-year  variation  may  be  mismodeled 
in  the  analytical  expressions.  The  frequency  of  the  analytical  term  does  not  quite  match  that  seen 
in  the  plot,  and  there  is  also  a  considerable  difference  in  the  amplitudes  of  the  numerical  and  an¬ 
alytical  terms.  Additionally,  a  36-year,  long-periodic  effect  is  suggested  by  the  plot;  this  term  is 
not  in  the  current  analytical  model.  Importantly,  however,  the  bound  on  eccentricity  variations,  as 
determined  by  the  analytical  model,  is  still  larger  than  observed  variations.  Consequently,  the  lunar 
model,  while  functioning  well  enough  for  present  purposes,  can  certainly  bear  further  investigation. 
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Figure  4.17:  Lunar  Third  Body  Induced  Variations  in  Radial  Direction. 


4.4.3  Summary  of  Lunar  Third  Body  Perturbation  Effects 

The  lunar  third  body  effects  on  radial  excursion  from  the  analytical  model  are  plotted  in 
Figure  4.17.  Radial  excursion  appears  to  increase  linearly  with  altitude.  In  the  altitude  region  of 
interest  in  this  study,  the  net  result  is  about  40-42  km.  There  are  no  eccentricity  nor  inclination 
dependence  in  the  major  effects,  seen  in  Table  4.7.  Short-periodic  terms  are  numerous.  There  are 
dozens  of  terms  with  the  same  order  of  magnitude,  and  no  one  term  clearly  dominates.  The  ones 
listed  in  Table  4.7  have  larger  amplitudes  for  the  particular  case  chosen  for  validation.  There  are 
also  instances  in  which  multiple  terms  of  the  same  frequency  with  have  large  amplitudes  individually 
but  with  opposite  signs  such  that  they  cancel  out  one  another.  Again,  the  periods  listed  in  the 
table  corresponds  to  the  test  case.  Lunar  effects  are,  in  fact,  the  largest  perturbations  aside  from 
deep  resonance  near  GEO  and  from  cases  involving  large  eccentricity  or  area-to-mass  ratios. 
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Table  4.7:  Analytical  Expressions  for  Significant  Lunar  Third  Body  Terms. 


Term 

Period 

Semimajor  Axis 

SP 

4.1015  •  10 ~16a5/2  cos  i$  cos  D3  sin  SI3  cos[2(03  -  A)] 

12  hr 

2.0508  •  10~16a5/2  cos2  *3  cos  D3  sin  ^3  cos[2(03  -  A)] 

12  hr 

1.3337  •  10-22a7/2  [sinf23cos(03  -  A)  +  cos^3sin(03  —  A)] 

24  hr 

LP 

none 

Eccentricity 

SP 

3.0762  •  10~16a3/2  cos  cos  Q3  sin  ^3  sin(2#3  —  A) 

24  hr 

LP 

A&nzaezj (96ni?3)  [5(1  —  cos  *3)  -  (1  +  cos  *3)]  sin  0J3 

9.1  yr 

4.5  Solar  Radiation  Pressure 

The  SRP  disturbing  potential,  as  noted  in  Section  3.7,  is  nearly  identical  to  that  of  third 
body  perturbations.  It  is  no  surprise  then  that  the  methodology  for  SRP  analysis  is  very  similar 
that  for  solar  third  body.  Unlike  the  other  perturbation  models,  this  one  required  no  patchwork; 
all  significant  terms  were  captured  in  the  basic  model. 

4.5.1  SRP  Analysis 

First,  in  order  to  obtain  a  disturbing  potential  function  for  SRP,  the  assumption  that  the 
satellite  is  always  sunlit  was  necessary.  For  GEO  altitudes,  this  is  not  unreasonable  and  will 
be  justified  in  the  next  section.  The  truncation  decisions  for  the  infinite  series  as  described  in 
Section  4.3.1  apply  here  without  changes,  though  the  decision  for  the  maximum  parallax  factor 
employed  is  no  longer  based  on  truncation  algorithms  embedded  in  DSST.  The  reason  for  this  is 
that  DSST  uses  a  shadow  model  similar  to  the  Cowell  integrator  and  does  not  assume  that  objects 
are  always  sunlit.  Instead,  the  maximum  parallax  factor,  Af  =  2,  is  implemented  simply  from 
experience  with  the  solar  model. 
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The  perturbing  source  is  again  assumed  to  be  in  a  circular  orbit,  and  so  Rs  is  constant  and 
equal  to  the  mean  distance  between  Earth  and  the  Sun.  The  factors  comprising  the  parameter  T, 
as  given  in  Equation  3.75,  are  fully  adjustable  in  the  SRP  model.  Default  values  for  validation  and 
application  purposes  are 

CR  =  1.2 

A/m  =  1  •  10-8  km2/kg 

C/c  =  4.51  •  10-3  kg/(s2km) 

Rs  =  1AU 

The  first  order  eccentricity  and  zero  satellite  inclination  assumptions  in  the  variation  equations 
met  with  no  adverse  consequences.  Similarly,  the  direction  cosines,  Equation  4.9,  with  the  values 
assigned  in  Equations  4.10  and  4.11,  worked  well  for  this  model.  Finally,  analytical  integration 
followed  the  same  lines  as  the  solar  third  body,  and  the  validation  process  shows  that  all  significant 
effects  are  reflected  in  this  SRP  model. 

4.5.2  SRP  Model  Validation 

The  validation  case  is  the  same  one  used  for  the  other  three  perturbation  models.  The 
significant  effects  due  to  SRP  are  clear-cut,  and  they  are  listed  in  Table  4.8.  Correlation  between 
analytical  and  numerical  results  are  very  good.  Short- periodic  effects  are  small  in  semimajor 
axis  and  negligible  in  eccentricity,  but  agreement  between  the  analytical  and  numerical  model  is 
excellent.  Plots  for  SP  variations  are  unnecessary  because  of  their  minimal  contribution  to  overall 
SRP  effects. 

Figure  4.18  shows  that  long-periodic  variations  in  semimajor  axis  are  on  the  submeter  level. 
The  twice-yearly  blips  in  the  plot  are  due  to  eclipsing  effects.  This  small  variation  in  semimajor  axis 
validates  the  always-sunlit  assumption  that  enables  the  formulation  of  a  SRP  disturbing  potential 
function.  The  plot  was  generated  by  DSST  AOG  using  a  SRP  with  model,  as  was  the  plot  in 
Figure  4.19  for  the  LP  variations  in  eccentricity.  A  simple  1-year  effect  in  eccentricity  is  captured 
in  this  plot.  This  term  represents  the  only  large  SRP  effect. 
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Table  4.8:  Significant  Periodic  Effects  from  SRP  Perturbations. 


Numerical  Amplitude  Analytical  Amplitude  Period 


Semimajor  Axis 


SP  20  m  20  m  24  hr 


LP  none 


Eccentricity 

SP  le-08  le-08  12  hr 


LP  1.165  e-04  1.3  e-04  1  yr 


Figure  4.18:  SRP  Long-Periodic  Variations  in  Semimajor  Axis. 


Time  (days) 


Figure  4.19:  SRP  Long  Periodic  Variations  in  Eccentricity. 


11.12 


Figure  4.20:  SRP  Induced  Variations  in  Radial  Direction. 

4.5.3  Summary  of  SRP  Perturbation  Effects 

SRP  perturbation  effects  in  the  radial  direction  are  plotted  in  Figure  4.20.  As  noted  before, 
the  short-periodic  variations  are  small,  in  fact  accounting  for  only  about  40  m  out  of  total  radial 
excursions  on  the  order  of  11  km.  The  rather  interesting  curves  in  the  plot  are  due  to  manual 
truncation  of  the  total  variations  in  semimajor  axis  and  eccentricity,  as  given  by  the  analytical 
model. 

Analytical  expressions  for  the  significant  terms  in  the  SRP  model  are  included  in  Table  4.9. 
The  argument  of  latitude  0©  of  the  perturbing  body  is  the  same  as  0$  in  Section  4.3.1.  Likewise, 
its  secular  rate  is  calculated  by  Equation  4.13.  The  only  large  effect,  the  yearly  long  periodic 
in  eccentricity,  dominates  the  SRP  model.  The  magnitude  of  this  term  is  linear  with  the  factor 
CrA/tu.  In  Figure  4.20,  CrA/tu  =  1.2  •  10-8  km2/kg  for  which  the  radial  excursion  is  roughly 
11  km.  For  cases  in  which  this  factor  is  large,  SRP  could  be  the  most  important  perturbation  effect. 
None  of  the  significant  terms  display  dependence  on  eccentricity  nor  inclination. 
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Table  4.9:  Analytical  Expressions  for  Significant  SRP  Terms. 


Term 

Period 

Semimajor  Axis 

SP 

1.3571  •  10-19a3/2T sin  (A  -  0©) 

24  hr 

5.8402  •  10-21<z3/2T sin  (A  +  0©) 

24  hr 

LP 

none 

Eccentricity 

SP 

-3.3924  •  10-2V/2rcos(2A  -  0©) 

12  hr 

-1.4601  •  10-21a1/2T  cos(2A  +  0©) 

12  hr 

LP 

1.0616-  10_1V/2r  cos  0© 

1  yr 

4.6  Cumulative  Effects 

After  studying  the  individual  perturbation  models,  the  combined  effects  of  all  perturbation 
models  were  investigated.  The  cumulative  effects  are  based  on  total  variations  from  each  component 
which  are  then  summed.  Figure  4.21  summarizes  the  cumulative  effects  on  radial  excursion  over 
a  range  of  altitudes  above  GEO.  This  plot  uses  the  parameters  of  the  test  case  with  zero  satellite 
inclination.  It  is  seen  that  geopotential  perturbations  dominate  until  about  GEO+20  km,  at  which 
point  the  lunar  effects  become  the  principal  player.  By  GEO+lOO  km,  geopotential  effects  are  the 
smallest,  falling  below  the  SRP  curve.  The  cumulative  curve  seems  to  follow  the  general  pattern  of 
the  geopotential  curve  until  about  GEO+400  km.  While  geopotential  effects  continue  to  decrease, 
the  other  perturbation  effects  are  increasing  enough  to  allow  the  total  curve  to  begin  an  upward 
slope. 

Numerical  integration  results  for  the  validation  case  with  all  perturbation  models  turned 
on  are  plotted  in  Figures  4.22  and  4.23  for  zero  inclination  and  in  Figures  4.24  and  4.25  for  15° 
inclination  over  20  years.  The  total  semimajor  axis  variation  for  both  inclination  cases  is  about 
23  km.  The  dominant  effects,  at  a  =  GEO +  30  km,  are  the  resonance  terms.  For  eccentricity,  the 
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Altitude  above  GEO  (km) 


Total 

- Geopotential 

. -  Lunar 

- Solar 

SRP 


Figure  4.21:  Cumulative  Effect:  Variations  in  Radial  Direction. 


Figure  4.22:  Total  Perturbation  LP  Effects  in  Semimajor  Axis  at  Zero  Inclination. 
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Eccentricity 


Figure  4.23:  Total  Perturbation  LP  Effects  in  Eccentricity  at  Zero  Inclination. 


8 


Time  (days) 


Figure  4.24:  Total  Perturbation  LP  Effects  in  Semimajor  Axis  at  15  Degrees  Inclination 


Eccentricity 


8.40E-03  - 

8.20E-03 

8.00E-03 

7.80E-03 

7.60E-03 

7.40E-03 

7.20E-03 

7.00E*03 

Figure  4. 


Table  4.10:  Comparison  of  Models  at  Different  Altitudes. 


Altitude 

Numerical 

Analytical 

GEO 

104  km 

136  km 

GEO+30  km 

78  km 

95  km 

GEO+300  km 

67  km 

75  km 

zero  inclination  case  shows  total  variations  on  the  order  of  1.4  •  10-3,  and  1.55  •  10-3  for  the  15° 
case,  both  primarily  due  to  lunar  perturbations.  These  variations  result  in  total  radial  excursions 
of  no  more  than  78  km  for  the  zero  inclination  case  and  88  km  for  the  other  one.  The  numbers 
for  radial  excursion  are  estimated  from  these  DSST  AOG  plots  in  conjunction  with  Cowell  plots 
showing  SP  variations.  Total  variations  are  differences  between  the  largest  and  smallest  values  on 
the  plots.  The  analytical  model  as  plotted  in  Figure  4.21  estimates  the  total  value  at  90  km  for 
zero  satellite  inclination. 

Finally,  Table  4.10  summarizes  the  comparison  of  the  analytical  against  the  numerical 
models  at  three  different  altitudes  with  all  four  perturbation  sources  modeled.  The  quantities  are 
maximum  radial  excursion  for  cases  with  inclinations  as  large  as  15°.  As  anticipated,  the  analytical 
model  consistently  estimates  more  conservatively  than  numerical  results.  A  positive  trend  is  that 
agreement  between  the  models  improves  as  the  altitude  increases.  The  reason  for  this  is  unclear.  It 
has  been  seen  that  lunar  effects  dominate  at  higher  altitudes,  but  this  study  notes  that  the  analytical 
model  for  lunar  third  body  perturbations  contains  the  greatest  uncertainty.  This  seeming  paradox 
is  left  to  future  investigations. 
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Chapter  5 


Application  to  the  Reorbiting 
Problem 


The  study  of  collisions  in  SSO  and  its  attending  implications  in  the  determination  of  the 
reorbiting  distance  from  GEO  requires  two  tools.  One  is  a  breakup  model  which  is  described  in  the 
next  section.  The  other,  a  perturbation  analysis  model,  is  summarized  in  Chapters  3  and  4.  The 
application  of  these  tools  toward  the  reorbiting  problem  is  final  objective  of  the  research. 

5.1  Breakup  Model 

At  the  time  of  this  research,  a  validated  breakup  model  suitable  for  the  GEO  regime  is  not 
available.  However,  the  purpose  of  the  research  is  not  to  build  a  low- velocity  collision  model  but 
to  study  short-  and  long-term  orbital  parameters  in  conjunction  with  the  SSO  collisions  scenario. 
Hence,  the  only  alternative  is  to  find  a  reasonable  and  available  breakup  model.  Such  a  model  need 
not  be  specifically  tailored  for  nonhypervelocity  collisions,  but  the  model  needs  to  provide  values  for 
the  delta-velocity  (Au)  imparted  to  the  fragmentation  pieces.  A  value,  associated  with  an  average 
piece,  is  consequently  used  to  determine  the  new  orbit  of  that  piece.  Then  the  perturbation  model 
is  applied,  and  the  region,  in  terms  of  radial  distance  that  the  representative  piece  may  occupy,  is 
established.  The  reorbiting  problem  seeks  to  control  the  placement  of  that  region. 


The  first  candidate  examined  was  the  statistical  breakup  model  in  ODESI,  which  is  de¬ 
scribed  in  Section  2.3.1.  The  environmental  model  ODESI  employed  the  assumption  that  the 
imparted  Av  to  the  collision  fragments  is  equal  to  the  collision  velocity,  i.e.,  the  relative  velocity 
between  the  parent  objects  [43],  and  this  assumption  seemed  to  be  a  reasonable  starting  point.  It 
was  understood  at  the  outset  that  this  procedure  would  produce  inflated  values  for  Av.  However, 
after  the  numbers  were  calculated,  it  was  clear  that  this  assumption  is  unacceptable.  For  instance, 
a  creditable  worst  case  scenario  is  the  collision  between  one  object  at  near  zero  inclination  and  the 
other  at  15°.  (Cases  of  collisions  between  objects  with  ascending  node  180°  apart  are  not  treated  in 
this  study  because  they  represent  implausible  events.)  If  both  objects  are  assumed  to  be  in  circular 
orbits,  the  relative  velocity  is 

vrel  =  0.26105^  (5.1) 

where  (jl  is  the  gravitational  parameter,  and  a,  the  semimajor  axis,  is  allowed  to  vary.  Assuming 
then  that  the  Av  is  also  given  by  Equation  5.1  and  also  that  it  is  applied  in  plane  and  at  a  full  180° 
from  the  velocity  direction  of  one  of  the  parent  objects,  the  new  semimajor  axis  and  eccentricity 
of  the  fragment  can  be  determined.  These  orbital  elements  in  turn  allow  for  the  calculation  of 
the  minimal  radial  distance  of  the  fragment  before  perturbations  affects  its  orbit.  Under  these 
assumptions,  a  collision  occurring  at  70,000  km  above  GEO,  produces  breakup  pieces  that  cross 
GEO  altitudes.  These  results  are  unrealistic  due  to  the  large  value  for  imparted  Av. 

Another  possibility  involves  somehow  scaling  the  results  from  a  hypervelocity  collision 
model.  This  may  be  an  abuse  of  the  model  as  it  was  not  intended  for  low-speed  encounters, 
but  on  the  flip  side,  research  has  not  conclusively  shown  that  such  models  are  indeed  inappropriate 
out  of  the  hypervelocity  realm.  This  alternative  approach  begins  with  the  data  generated  from  the 
simulation  of  a  collision  between  a  satellite  in  GEO  and  a  debris  piece  in  geosynchronous  transfer 
orbit  [21].  The  relative  velocity  is  on  the  order  of  1.4  km/s.  This  simulation  is  based  on  a  hyper- 
velocity  model,  and  conservation  of  mass  and  energy  is  confirmed  [17].  Subjective  approximations 
from  graphed  data  of  the  spatial  distribution  of  collision  debris  provides  orbital  information  for  a 
selected  average  piece,  which  is  discussed  in  further  detail  later.  The  semimajor  axis  is  estimated, 
and  with  the  assumption  of  a  full  retrograde,  in-plane  impulse,  the  orbit  is  taken  as  being  elliptical 
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with  the  apogee  height  at  the  GEO  parent  object’s  initial  altitude.  The  difference  between  the 
velocity  of  the  parent  object’s  circular  orbit  and  the  velocity  of  the  fragment  at  apogee  is  the  Av 
imparted  by  the  collision.  This  is  calculated  to  be  20  m/s.  A  linear  scaling  function  is  employed 
to  adjust  the  Av  for  different  impact  velocities: 

Av  =  (5-2) 

=  0.014286vre;  (5.3) 

All  quantities  are  in  units  of  km/s. 

The  average  piece  chosen  for  this  approach  is  the  median  data  point  in  the  set  of  debris 
pieces  whose  semimajor  axis  is  smaller  than  that  of  the  GEO  parent  object.  Since  the  direction 
component  of  imparted  Av  is  completely  random,  roughly  half  of  the  fragmentation  debris  are  flung 
in  a  prograde  direction  while  the  rest  experience  retrograde  impulses.  The  debris  size  distribution 
of  these  two  sets  should  be  nearly  identical.  The  relation  between  debris  size  and  imparted  Av  is 
more  complicated.  Naturally,  the  smaller  pieces  are  more  likely  to  experience  larger  Av.  Therefore, 
a  reasonable  claim  may  be  made  that,  as  the  chosen  piece  is  a  median  point,  roughly  one  quarter  of 
the  collision  fragments  will  be  orbits  lower  than  the  representative  piece,  regardless  of  the  breakup 
mass  distribution,  and  that  these  fragments  will  be  the  smaller  fragments. 

With  the  relation  in  Equation  5.3,  the  assumption  on  the  direction  of  the  impulse,  and 
two-body  orbital  equations,  the  orbit  geometry  of  the  fragment  can  be  calculated  for  given  collision 
altitudes.  For  collisions  between  objects  in  circular  orbits  of  the  same  semimajor  axis  and  separated 
by  15°  in  inclination,  data  for  varying  collision  altitudes  are  presented  in  Table  5.1.  Eccentricity  of 
the  fragment  in  its  initial  orbit  after  the  breakup  is  0.00744  for  all  reorbiting  altitudes,  and  r 
is  the  minimal  radial  distance,  based  on  two-body  perigee  distance.  For  this  study,  the  collision 
altitude  is  taken  to  be  equal  to  reorbiting  altitudes.  A  more  moderate  case  is  summarized  in 
Table  5.2  for  collisions  between  circular  orbits  separated  by  7.5°.  The  eccentricity  of  the  fragment 
is  smaller,  0.00373,  due  to  the  smaller  imparted  Au.  In  both  cases,  semimajor  axis  for  GEO  is 
42,164  km. 

It  is  important  to  note  that  the  minimal  radial  distances  in  Tables  5.1  and  5.2  are  based  only 
on  two-body  equations  without  perturbations  effects  taken  into  account.  The  numbers  themselves, 
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Table  5.1:  Data  From  Collisions  Between  Objects  15  Degrees  Apart. 


Reorbiting  Altitude 

Av 

New  Orbit  of  Fragment 

[km] 

[km/s] 

a  [km] 

Train  [km] 

300 

41836 

400 

41935 

500 

0.0114 

42349 

42033 

600 

42448 

42132 

633 

42481 

42164 

700 

42230 

800 

0.0114 

42647 

42329 

900 

0.0113 

42746 

42428 

Table  5.2:  Data  From  Collisions  Between  Objects  7.5  Degrees  Apart. 


Reorbiting  Alt 

Av 

New  Orbit  of  Fragment 

[km] 

[km/s] 

a  [km] 

r min  [km] 

300 

0.00573 

42306 

42148 

317 

0.00572 

42323 

42165 

400 

0.00572 

42406 

42247 

500 

0.00571 

42505 

42347 

600 

0.00571 

42605 

42446 

700 

0.00570 

42705 

42545 
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however,  are  consistent  with  expectations.  It  is  no  surprise  that  collisions  occurring  at  300  km 
above  GEO  will  produce  debris  that  crosses  GEO.  This  breakup  model  cannot  be  validated  by  any 
available  means  at  this  time  due  to  the  lack  of  empirical  data.  However,  a  methodology  is  established 
such  that  a  different  model  may  be  substituted  for  this  one,  should  it  prove  to  be  better.  This  tool 
can  now  be  applied  in  conjunction  with  the  perturbation  models  to  study  reorbiting  altitudes. 


5.2  Reorbiting  Altitudes 

Six  cases  were  studied  in  which  the  definition  of  the  GEO  region  and  the  inclination  between 
the  colliding  objects  varied.  The  GEO  region  is  defined  in  three  ways: 

•  Case  1:  Restricted  to  GEO  altitude,  42,164  km 

•  Case  2:  Covers  the  region  traversed  by  uncontrolled  GEO  objects,  GEO+104  km 

•  Case  3:  Includes  drifting  and  operational  regions  and  a  modest  buffer,  GEO+200  km 

The  first  case  represents  the  absolutely  minimal  reorbiting  requirement:  keeping  most  debris  clear 
of  operational  satellites.  In  reality,  functional  GEO  satellites  occupy  various  altitudes  around  the 
GEO  altitude,  but  this  case  is  meant  to  illustrate  an  extreme  definition  of  the  GEO  region.  The 
next  case  involves  a  moderate  definition.  The  region  occupied  by  drifting  GEO  satellites  was 
determined  by  long-term  semianalytical  integration,  modeling  a  4x4  JGM-2  gravity  field,  lunisolar 
perturbations,  and  SRP  effects  for  A/m  —  0.01  m2/kg.  The  last  definition  for  the  GEO  region  is 
conservative  by  some  standards  because  it  includes  a  buffer,  the  value  for  which  is  generally  not 
based  upon  physics  or  facts  but  is  set  by  individual  estimations. 

For  each  of  the  three  GEO  region  cases,  two  collision  scenarios  are  examined.  Since  uncon¬ 
trolled  objects  in  the  GEO  regime  generally  have  inclinations  between  0°  and  15°,  a  worst  case  and 
a  more  typical  case  are  studied.  The  worst  case  is  defined  by  collisions  between  objects  separated 
by  the  full  15°.  Recognizing  that  this  is  an  extreme  case  and  that  preventive  measures  might  be 
geared  more  toward  typical  situations,  collisions  between  objects  separated  by  7.5°  in  inclination 
are  also  surveyed.  There  are  then  six  distinct  cases  discussed  in  the  following  sections. 
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Table  5.3:  GEO  Region  Case  1  with  Collisions  Between  Objects  15  Degrees  Apart. 


GEO  region 

42164  km 

Before  collision 

Semimajor  axis 

42884  km 

Reorbiting  distance 

720  km 

After  collision 


Semimajor  axis 

42567  km 

Eccentricity 

0.0074447 

Perturbations 

A  a  [km] 

Ae 

Geopotential 

1.9715 

0.0001286 

Solar  3rd  body 

1.0043 

0.0004434 

Lunar  3rd  body 

2.4114 

0.0009382 

SRP 

0.0506 

0.0003199 

Total 

5.4378 

0.0018601 

Min  radial  distance 

42167  km 

5.2.1  15  Degrees  Inclination  Case 

Each  definition  of  the  GEO  region  specifies  an  altitude  that  the  representative  fragment 
should  not  dip  below.  Using  the  breakup  model  described  in  Section  5.1  and  the  analytical  pertur¬ 
bation  models  developed  in  this  research,  the  collision  altitude,  which  is  also  the  reorbiting  altitude, 
is  adjusted  manually  until  the  models  show  that  the  fragment  will  not  enter  the  specified  region. 
A  5  km  tolerance  is  allowed.  Tables  5.3,  5.4,  and  5.5  summarize  the  results  for  the  three  differently 
defined  GEO  regions. 

The  semimajor  axis  before  the  collision  is  the  parameter  sought.  The  reorbiting  distance 
is  the  difference  of  this  value  and  42,164  km.  Recall  that  the  breakup  model  assumes  that  the 
parent  objects  are  initially  in  circular  orbits  before  the  collision.  The  perturbation  models  uses  the 
semimajor  axis  and  eccentricity  calculated  by  the  breakup  model  and  in  turn,  gives  the  maximum 
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Table  5.4:  GEO  Region  Case  2  with  Collisions  Between  Objects  15  Degrees  Apart. 


GEO  region 

42268  km 

Before  collision 

Semimajor  axis 

42989  km 

Reorbiting  distance 

825  km 

After  collision 


Semimajor  axis 

42671  km 

Eccentricity 

0.0074447 

Perturbations 

Aa  [km] 

Ae 

Geopotential 

1.6133 

0.0001278 

Solar  3rd  body 

1.0141 

0.0004468 

Lunar  3rd  body 

2.4358 

0.0009429 

SRP 

0.0510 

0.0003203 

Total 

5.1143 

0.0018379 

Min  radial  distance 

42270  km 
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Table  5.5:  GEO  Region  Case  3  with  Collisions  Between  Objects  15  Degrees  Apart. 


GEO  region 

42364  km 

Before  collision 

Semimajor  axis 

..  . 

43089  km 

Reorbiting  distance 

925  km 

After  collision 


Semimajor  axis 

42771  km 

Eccentricity 

0.0074447 

Perturbations 

A  a  [km] 

Ae 

Geopotential 

1.3866 

0.0001272 

Solar  3rd  body 

1.0236 

0.0004501 

Lunar  3rd  body 

2.4593 

0.0009473 

SRP 

0.0513 

0.0003207 

Total 

4.9208 

0.0018453 

Min  radial  distance 

42368  km 
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variations  in  semimajor  axis,  A  a,  and  in  eccentricity,  Ae.  The  minimal  radial  distance  is  calculated 

by 

Train  =  ifineui  ~  Aft)  (1  —  ( 6new  +  Ae))  (S.4) 

where 

anew  =  the  fragment’s  semimajor  axis  after  the  collision 
enew  =  the  fragment’s  eccentricity  after  the  collision 

This  expression  is  based  upon  the  two-body  equation  for  periapse. 

Tables  5.3,  5.4,  and  5.5  indicate  that  the  reorbiting  distance  is  largely  driven  by  the  Av 
from  the  breakup  model.  Unfortunately,  there  is  much  uncertainty  in  this  model,  and  therefore,  the 
estimates  for  reorbiting  distance  should  be  treated  carefully.  Of  the  perturbation  effects,  several 
trends  are  noted.  For  all  disturbing  forces  except  geopotential,  the  associated  Aft  and  Ae  increase 
with  altitude.  The  geopotential  effects  decrease  with  increasing  altitude,  as  expected  from  the 
results  in  Chapter  4.  Increases  in  reorbiting  altitude  seem  to  hold  a  linear  relationship  with  the 
varying  definitions  for  the  GEO  protected  region.  If  some  faith  is  placed  in  the  breakup  model, 
then  to  protect  the  barest  GEO  region,  Case  1,  current  reorbiting  guidelines  need  to  be  doubled  at 
the  least.  The  economic  penalties  may  not  be  reasonable. 

5.2.2  7.5  Degrees  Inclination  Case 

The  same  procedure  is  repeated  for  the  more  moderate  collision  case  of  7.5°  inclination 
difference  between  the  objects.  Results  are  given  in  Tables  5.6,  5.7,  and  5.8.  As  expected,  the 
reorbiting  distances  are  smaller  in  these  cases.  However,  the  least  expensive  option,  Case  1,  still 
requires  a  reorbiting  distance  that  is  30%  larger  the  current  guidelines  of  300  km.  Previously 
noted  trends  in  the  models  are  applicable  here  as  well.  Unlike  in  the  extreme  cases  summarized 
in  Section  5.2.1,  the  reorbiting  requirements  shown  here  may  not  be  considered  unreasonable.  The 
most  conservative  case,  for  a  GEO+200  km  protected  region,  indicates  that  590  km  may  be  a 
sufficient  reorbiting  distance.  AUSSAT  at  one  time  suggested  reorbiting  up  to  1110  km  above 
GEO  [13]. 
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Table  5.6:  GEO  Region  Case  1  with  Collisions  Between  Objects  7.5  Degrees  Apart. 


GEO  region 

42164  km 

Before  collision 

Semimajor  axis 

42554  km 

Reorbiting  distance 

390  km 

After  collision 


Semimajor  axis 

42396  km 

Eccentricity 

0.0037338 

Perturbations 

A  a  [km] 

Ae 

Geopotential 

3.2381 

0.0001293 

Solar  3rd  body 

0.9773 

0.0002411 

Lunar  3rd  body 

2.3394 

0.0008070 

SRP 

0.0499 

0.0003193 

Total 

6.6046 

0.0014966 

Min  radial  distance 

42167  km 
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Table  5.7:  GEO  Region  Case  2  with  Collisions  Between  Objects  7.5  Degrees  Apart. 


GEO  region 

42268  km 

Before  collision 

Semimajor  axis 

42659  km 

Reorbiting  distance 

495  km 

After  collision 


Semimajor  axis 

42500  km 

Eccentricity 

0.0037338 

Perturbations 

A  a  [km] 

Ae 

Geopotential 

2.2938 

0.0001285 

Solar  3rd  body 

0.9870 

0.0002430 

Lunar  3rd  body 

2.3632 

0.0008115 

SRP 

0.0503 

0.0003197 

Total 

5.6943 

0.0015027 

Min  radial  distance 

42272  km 
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Table  5.8:  GEO  Region  Case  3  with  Collisions  Between  Objects  7.5  Degrees  Apart. 


GEO  region 

42364  km 

Before  collision 

Semimajor  axis 

42754  km 

Reorbiting  distance 

590  km 

After  collision 


Semimajor  axis 

42595  km 

Eccentricity 

0.0037338 

Perturbations 

Aa  [km] 

Ae 

Geopotential 

1.8349 

0.0001278 

Solar  3rd  body 

0.9958 

0.0002447 

Lunar  3rd  body 

2.3850 

0.0008157 

SRP 

0.0506 

0.0003200 

Total 

5.2663 

0.001508 

Min  radial  distance 

42366  km 

93 


Chapter  6 


Conclusions  and  Future  Work 


The  goal  of  this  research  appeared  clear-cut  at  the  outset  and  likewise,  the  path  to  achieve 
it  seemed  simple  at  first  glance.  At  the  beginning,  the  target  of  these  studies  looked  to  be  a 
single  number,  the  minimum  safe  reorbiting  distance.  Such  was  the  naivete  that  optimistically 
colored  early  decisions  and  projections.  The  low-velocity  breakup  models  available  at  the  time  that 
research  began  seemed  unsuitable,  but  there  was  confidence  that  better  models  would  be  ready 
in  the  near  future.  Two  years  later,  an  unvalidated  model  is  incorporated  into  this  study.  The 
reorbiting  problem,  as  addressed  in  Chapter  5,  did  not  have  a  single  number  as  its  answer.  Instead, 
summaries  filling  six  tables  attempt  to  describe  different  aspects  of  a  situation  fraught  with  many 
alternatives  and  uncertainties.  This  study  provides  an  initial,  and  perhaps  too  simplistic,  look  at 
the  effects  of  collisions  in  SSO.  Refinement  possibilities  are  numerous. 

6.1  Conclusions 

This  research  may  be  divided  into  three  components,  and  each  are  now  addressed  sepa¬ 
rately.  The  first  part  concerns  the  perturbation  models.  Though  analytical  methods  sometime 
produce  unwieldy  expressions,  the  gains  include  physical  insight,  and  specific  information  on  the 
relationship  between  individual  effects  and  orbital  parameters,  such  as  semimajor  axis,  eccentricity, 
and  inclination,  and  that  between  effects  and  physical  characteristics  of  the  satellite.  Furthermore, 
the  variation  expressions  are  applicable  to  a  whole  class  of  orbits.  All  of  this  is  unavailable  from 
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numerical  methods,  except  through  inference.  An  additional  advantage  comes  from  the  ease  with 
which  the  analytical  perturbation  model  interfaces  with  the  breakup  model.  Results  for  the  re- 
orbiting  problem  equaling  years  of  propagation  are  produced  in  mere  seconds.  Amplitudes  and 
frequencies  of  individual  effects  are  more  accurately  estimated  analytically  than  through  orbital 
element  history  plots  showing  multiple  effects  enmeshed  together. 

The  results  of  the  model  validation  demonstrate  that  the  majority  of  the  perturbation 
effects  are  accurately  captured  by  the  variation  equations.  Short-periodic  terms  are  generally 
modeled  better  in  the  analytical  models  since  they  are  better  suited  to  the  analytical  integration 
method  of  holding  certain  slower-varying  quantities  constant.  On  the  other  hand,  some  long-term 
propagations  indicate  that  further  investigation  into  the  long-periodic  effects  is  warranted.  In  spite 
of  this,  the  results  of  the  analytical  models  presented  in  this  study  do  bound  all  variations  over 
20-year  propagations.  In  fact,  Table  4.10  shows  that  the  analytical  models  do  their  best  in  the 
currently  recommended  SSO  region  in  the  vicinity  of  GEO+300  km. 

The  second  component  is  the  breakup  model.  From  the  earlier  descriptions  of  the  model 
employed,  it  is  easy  to  conclude  that  much  of  the  uncertainty  in  the  reorbiting  analysis  stems 
from  use  of  this  unvalidated  model.  While  it  produced  seemingly  reasonable  numbers,  estimating 
the  error  bounds  is  impossible  because  there  is  no  truth  model  for  comparison.  Nonhypervelocity 
breakup  dynamics  is  a  largely  unexplored  area;  the  first  steps  sometimes  must  be  taken  on  faith. 
A  better  model,  when  it  becomes  available,  can  easily  be  substituted  in  for  the  current  one. 

With  the  methodology  established  and  given  reliable  tools,  the  reorbiting  analysis  can  be 
accomplished  swiftly.  Table  6.1  summarizes  the  reorbiting  distance  for  the  six  cases  examined  in 
Chapter  5.  The  accuracy  of  this  analysis  hinges  greatly  on  the  breakup  model’s  accuracy.  The 
reorbiting  distance  is  driven  mostly  by  the  collision,  with  perturbations  adding  less  than  100  km. 
Since  definitions  for  the  desired  protected  GEO  region  may  vary,  the  relationship  between  the 
GEO  region  and  the  reorbiting  distance  is  explored.  It  appears  to  be  no  worse  than  a  linear 
relation.  The  importance  of  the  breakup  model  indicates  that  future  work  on  the  reorbiting  problem 
should  focus  there.  If  the  results  of  the  current  analysis  are  taken  at  face  value,  current  reorbiting 
practices  are  probably  not  sufficient  to  protect  the  GEO  region,  any  definition  of  it,  in  the  event 
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Table  6.1:  Summary  of  Reorbiting  Analysis. 


Collision  Case 

GEO  Region 

Reorbiting  Distance 

15  Degrees 

GEO  only 

720  km 

15  Degrees 

GEO+104  km 

825  km 

15  Degrees 

GEO+200  km 

925  km 

7.5  Degrees 

GEO  only 

390  km 

7.5  Degrees 

GEO+104  km 

495  km 

7.5  Degrees 

GEO+200  km 

590  km 

of  collisions.  Furthermore,  the  safe  reorbiting  distances,  as  calculated  in  this  study,  are  probably 
beyond  reasonable  cost  penalties  in  the  15°  cases.  The  7.5°  cases  may  also  be  considered  extreme. 
Current  technology  may  not  permit  consideration  of  collisions  in  SSO. 

6.2  Future  Work 

If  the  analytical  approach  is  continued,  improvements  can  be  made  in  three  particular  areas. 
The  first  concerns  the  geopotential  and  lunar  third  body  models.  As  stated  previously,  long-periodic 
effects  on  eccentricity  are  suspected  to  be  missing  or  mismodeled  in  the  current  analytical  model. 
Based  on  eccentricity  histories  from  semianalytical  orbit  propagations  for  the  lunar  case,  effects 
with  very  long  periods,  about  36  years,  are  present  but  not  included  in  the  current  model.  The 
long-periodic  term  in  the, geopotential  model  has  notably  poor  correlation  between  analytical  and 
semianalytical  estimates  of  the  amplitude  and  frequency.  Additionally,  nonzero  inclination  terms 
are  generally  ignored,  except  for  a  few  of  the  largest  terms.  Recovery  of  these  terms,  perhaps 
especially  in  the  solar  model,  would  help  to  refine  the  perturbation  models.  Finally,  the  singularity 
of  the  geopotential  model  near  GEO  altitudes  should  be  addressed,  if  only  for  completeness.  The 
handling  of  resonance  terms  which  are  responsible  for  this  difficulty  require  a  fresh  approach. 
However,  the  primary  issue  of  this  research  concerns  the  SSO  region  so  future  investigations  need 
not  focus  too  heavily  on  perturbations  of  GEO  which  has  also  been  studied  by  many. 
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A  completely  numerical  approach  to  perturbation  analysis  would  require  longer  propagation 
capabilities  in  order  to  view  longer  periodic  effects.  To  shorten  integration  times,  double  averaging 
theories  could  provide  a  great  advantage  [9].  A  numerical  study  of  the  reorbiting  problem  also 
requires  a  search  algorithm  to  confirm  maximum  variation  bounds. 

There  is  much  room  for  improvement  in  the  low-speed  collision  model.  A  better  model 
would  provide  typical  fragment  characteristics,  such  as  average  area-to-mass  ratios,  as  well  as  more 
accurate  mass  and  velocity  distributions.  For  validation  of  the  model,  empirical  data  from  ground 
tests  is  essential.  On-orbit  breakup  data  would  require  more  sensitive  detection  capabilities  than 
the  current  1-meter  diameter  limit.  Due  to  the  great  importance  of  the  breakup  model  in  the 
reorbiting  analysis,  emphasis  in  future  work  is  placed  in  this  area. 

Finally,  the  reorbiting  analysis  could  be  improved  by  incorporating  the  breakup  and  pertur¬ 
bation  models  to  form  one  easy-to-use  analysis  tool.  This  in  turn  can  be  a  part  of  a  comprehensive 
GEO  environmental  model.  The  user  would  be  able  to  set  his  own  definition  of  the  GEO  region, 
select  any  combination  of  perturbation  models,  and  even  determine  the  collision  scenario,  if  infor¬ 
mation  on  the  orbits  of  objects  predicted  for  a  near-approach  is  available.  Though  the  analytical 
perturbation  models  address  general  classes  of  orbits,  the  breakup  model  should  be  able  to  evaluate 
specific  cases  and  therefore,  be  of  use  in  risk  assessments  as  well  as  debris  mitigation  studies. 
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